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Counterions at charged rod-like polymers exhibit a condensation transition at a critical tempera- 
ture (or equivalently, at a critical linear charge density for polymers), which dramatically influences 
various static and dynamic properties of charged polymers. We address the critical and universal 
aspects of this transition for counterions at a single charged cylinder in both two and three spatial 
dimensions using numerical and analytical methods. By introducing a novel Monte-Carlo sampling 
method in logarithmic radial scale, we are able to numerically simulate the critical limit of infinite 
system size (corresponding to infinite-dilution limit) within tractable equilibration times. The criti- 
cal exponents are determined for the inverse moments of the counterionic density profile (which play 
the role of the order parameters and represent the inverse localization length of counterions) both 
within mean-field theory and within Monte-Carlo simulations. In three dimensions, we demonstrate 
that correlation effects (neglected within mean-field theory) lead to an excessive accumulation of 
counterions near the charged cylinder below the critical temperature (condensation phase), while 
surprisingly, the critical region exhibits universal critical exponents in accord with the mean-field 
theory. Also in contrast with the typical trend in bulk critical phenomena, where fiuctuations are 
strongly enhanced in lower dimensions, we demonstrate, using both numerical and analytical ap- 
proaches, that the mean-field theory becomes exact for the 2D counterion-cylinder system at all 
temperatures (Manning parameters), when number of counterions tends to infinity. For finite par- 
ticle number, however, the 2D problem displays a series of peculiar singular points (with diverging 
heat capacity), which refiect successive de- localization events of individual counterions from the 
central cylinder. In both 2D and 3D, the heat capacity shows a universal jump at the critical point, 
and the energy develops a pronounced peak. The asymptotic behavior of the energy peak location 
is used to locate the critical point, which is also found to be universal and in accordance with the 
mean-field prediction. 

PACS numbers: 64.60.Fr, 61.20.Ja, 82.35.Rs, 87.15.-v 



I. INTRODUCTION 

Electrostatics of charged polymers is often dominated 
by small oppositely charged ions (counterions), which 
maintain the global electroneutrality of charged solu- 
tions. Many charged polymers, such as tubulin, actin 
and DNA are stiff and may be represented by straight 
cylinders (on length scales smaller than the persistence 
length). Neglecting many- ion effects, a single counte- 
rion is attracted by an electrostatic potential that grows 
logarithmically with the radial distance from the cen- 
tral cylinder axis. But since the counterion confine- 
ment entropy also shows a logarithmic size dependence, 
it was suggested early by Onsager (Ij that a counte- 
rion delocalization transition occurs at a critical cylin- 
der charge or equivalently, at a critical temperature. 
Onsager 's argument, which is strictly valid for a sin- 
gle particle, was soon corroborated by mean-field studies 
[1 11 i, 11 IE IE Si, IS ES 111 El 111, which demonstrate 
that a charged cylinder can indeed bind or condense a fi- 
nite fraction of counterions below a critical temperature 
(and even in the limit of infinite system size with no con- 
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fining boundaries), while above the critical temperature, 
all counterions de-condense and diffuse to infinity. 

This counterion- condensation transition (CCT) dra- 
matically affects a whole number of static and dynamic 
quantities as observed in recent experiments on charged 

polymers 0,ll0,|ll|ll|ll|ll|ll|lSl23: upon con- 
densation, the bare polymer charge is screened leading, 
for instance, to a significant reduction in electrophoretic 
mobility 0, |23| and conductivity of polymers E3; it 
also triggers striking static properties such as counterion- 
induced attraction between like-charged polymers, which 
gives rise to compact phases of F-actin 21j and DNA [22| . 
Since its discovery, the CCT has been at the focus of nu- 
meric^lMMMM and analytica^l Illj2l^ M 
111 111 "3?,", ^,^111 111 lag, "iff, "ir^ stud- 
ies. Under particular dispute has been the connection be- 
tween CCT and the celebrated Kosterlitz-Thouless tran- 
sition of logarithmically interacting particles in two di- 
mensions 31, 45, 46, 47] . 

The CCT at charged cylinders is regulated by a dimen- 
sionless control parameter, ^ = q^BT, known as the Man- 
ning parameter J,], which depends on the linear charge 
density of the cylinder, —re, charge valency of counteri- 
ons, +q, and the Bjerrum length £b = e"^ / {AneeokBT) 
accommodating the ambient temperature T and the 
medium dielectric constant s. The Manning parame- 
ter plays the role of the inverse rescaled temperature 
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and can be varied experimentally by changing the lin- 
ear charge density (using synthetic chains or various pH) 
d 111 El 13 or by varyinff the dielectric media (mix- 
ing different solve nts) iVA l20l . According to mean- field 
theory ^1. 2. 9. 10. iTlL[l^ll3| . condensation occurs above 
the critical value = 1- In experiments, the critical 
Manning parameter appears to be about unity, but large 
deviations have also been reported f9, 20, 4SJ, and the 
precise location of the critical point is still debated . 

On the other hand, it is known that the critical tem- 
perature may in general be influenced by correlations and 
fluctuations, which are not captured within the mean- 
field theory ^9|. These effects typically cause devia- 
tions from mean-field predictions in both non-universal 
and universal quantities below the upper critical dimen- 
sion. Surprisingly, the mean-field prediction for the CCT 
threshold, ^c, has not been questioned in literature and 
apparently assumed to be exact. Likewise, the existence 
of universal scaling relations and critical (scaling) expo- 
nents associated with the CCT has not been addressed, 
neither on the mean-field level nor in the presence of cor- 
relations. 

Our chief goal in this paper is to address the follow- 
ing issues: i) what is the exact threshold of the CCT, 
^c, and ii) what are the critical exponents associated 
with this transition both in three and two spatial di- 
mensions. We shall also address the type of singularities 
that emerge in thermodynamic quantities as the CCT 
criticality sets in. To establish a systematic investiga- 
tion of the correlation effects, we employ Monte-Carlo 
simulations for counterions at a single charged cylinder 
using a novel sampling method (centrifugal sampling), 
which is realized by mapping the radial coordinate to a 
logarithmic scale. This enables us to investigate the crit- 
ical limit of infinite system size (that is when the outer 
boundaries confining counterions tend to infinity) within 
tractable equilibration times in the simulations. The im- 
portance of taking very large system sizes becomes evi- 
dent by noting that lateral finite-size effects, which mask 
the critical unbinding behavior of counterions, depend on 
the loqarithm of s yste m size in the cylindrical geometry 
lllllSIIillllllMlMIIllslllllMliilli^^ caus- 
ing a quite weak convergence to the critical infinite-size 
limit. 

Our simulations provide the first numerical results for 
the asymptotic critical behavior of CCT and system- 
atically incorporate correlation effects (a brief report 
of some of our results has been presented previously 
^3|)- The relevance of electrostatic correlations is in 
general identified by a dimensionless coupling parame- 
ter, S = 2TTq^i^as with = t/{2ttR) being the surface 
charge density and R the radius of the cylinder. The 
mean- field theory represents the limit S ^ [H^ 
while in the converse limit of strong coupling, S ^ 1, cor- 
relations become significant and typically lead to drastic 
changes [s^, 113, Hj, H^. In order to investigate scal- 
ing properties of the CCT in various regimes of the cou- 
pling parameter, we focus on the inverse moments of the 



counterionic density profile, which play the role of the 
"order parameters" for this transition and represent the 
mean inverse localization length of counterions. Using a 
combined finite-size-scaling analysis with respect to both 
lateral size of the system and the number of counterions, 
we show that the order parameters adopt scale-invariant 
forms in the vicinity of the critical point. The critical 
exponents associated with the reduced temperature and 
the size parameters are determined both within the simu- 
lations and also analytically within two limiting theories 
of mean field and strong coupling. As a main result, we 
find that the critical exponents of the CCT are univer- 
sal (that is independent of the coupling parameter varied 
over several decades 0.1 < S < 10^) and appear to be in 
close agreement with the mean-field prediction. Surpris- 
ingly, we find that the critical Manning parameter is also 
universal and given by the mean-field value = 1- The 
transition threshold, ^c, is determined with high accuracy 
from the asymptotic behavior of the location of a singu- 
lar peak that emerges in average internal energy of the 
system. The excess heat capacity is found to vanish at 
small Manning parameters (de-condensation phase) and 
exhibits a universal jump at the transition point indi- 
cating that the CCT may be regarded as a second-order 
phase transition as also suggested in a previous mean- 
field study |3^. 

As will be shown, the validity of mean-field predictions 
in 3D breaks down as the Manning parameter increases 
beyond the critical value (i.e. in the condensation phase), 
where inter-particle correlations become significant at 
large couplings. This leads to an enhanced accumulation 
of counterions near the cylinder surface and a crossover 
to the strong-coupling theory predictions 52, ,53^ 5(^ . 

In order to bring out possible role of fiuctuations, we 
also study the CCT in a 2D counterion-cylinder sys- 
tem (equivalent to a 3D system composed of a central 
charged cylinder and parallel cylindrical "counterions" 
with logarithmic Coulomb interactions, as may be appli- 
cable to an experimental system of oriented cationic and 
anionic polymers, e.g., DNA with polylysine jl^)- For 
finite number of counterions, a peculiar series of singular 
points emerge that refiect delocalization events of indi- 
vidual counterions as the Manning parameter varies. For 
increasing particle number, the singular points tend to 
merge and eventually in the thermodynamic limit, the 
2D results tend to universal values determined by the 
mean-field theory. Therefore, in contrast to the typical 
trend in critical phenomena, the CCT in 2D is found to 
be in exact agreement with the mean-field theory for the 
whole range of Manning parameters (or temperatures) 
when the number of counterions tends to infinity. As 
will be shown, the simulation results in 2D can be re- 
produced using an approximate analytical approach. A 
more systematic method has been developed recently [t^ 
supporting the present analytical results. 

The organization of the paper is as follows: in Sec- 
tions IIIIVII we focus on the counterion-cylinder system 
in three spatial dimensions. Our model is introduced in 
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Sectional where we shall also outline the general method 
proposed for the investigation of the CCT. In Section lTm 
we derive the scaling relations for order parameters and 
determine the asymptotic behavior of thermodynamic 
quantities within the mean-field theory (which is valid 
in all dimensions). In Section llVI analytical results are 
obtained within the strong-coupling theory. The numer- 
ical analysis of the CCT for various coupling strength 
will be presented in Sections IVII and IVIII in three and 
two dimensions, respectively. 



II. COUNTERION-CONDENSATION 
TRANSITION (CCT) IN THREE DIMENSIONS 

A. Cell model for charged rod-like polymers 

We consider a primitive cell model |^ |^ |^ , which 
consists of a single charged cylinder of radius R and 
point-like neutralizing counterions of charge valency ~\-q 
that are confined laterally in an outer (co-axial) cylin- 
drical box of radius D-see Figure ^ The cylinder has 
infinite length, L, and a uniform (surface) charge distri- 
bution, — tT(x)e, where cr(x) = as5{r — R). (Note that q 
and CTs are given in units of the elementary charge, e, and 
are positive by definition.) The cylinder is assumed to be 
rigid and impenetrable to counterions and the dielectric 
medium is represented by a uniform dielectric constant, 
£. In the three dimensions, electric charges interact via 
bare Coulombic interaction 



W3d(x) = l/|x|. 



(1) 



The electroneutrality condition holds globally inside 
the cell and entails the relation 



qN = tL, 



(2) 



where N is the number of counterions per cell and 
T = 2'KRfJs represents the linear charge density of the 
cylinder. The system is described by the Hamiltonian 

= g 2^ W3D(Xi - Xj ) 

N 

- g^B^ / W3D (x - Xi) cr(x) dx 

H- ^ y a(x)«3D(x-x')(T(x')dxdx', (3) 

which comprises mutual repulsions between counterions 
located at {x^} (first term), the counterion-cylinder at- 
traction (second term) and the self-energy of the cylinder 
(last term). It can be written as 



n 



N 



knT 



N 

= q^iB «3d(x. - X,) + 2^ ^ In (^) + Co, (4) 



- 2D 
\-2R- 




FIG. 1: The three-dimensional model consists of a charged 
cylinder of infinite length, H, and its neutralizing counterions 
confined in an outer cylindrical box (see the text for param- 
eters) . 



where ^ is the Manning parameter of the system P, Q , 



(5) 



with £b = e'^ / (AnseokBT) being the Bjerrum length 
(in water and at room temperature £-b — TA), and 
ri = (xf -t- yf y^^ being the radial coordinate of the i-th 
counterion from the cylinder axis, which coincides with 
z-axis. The additive term Co in Eq. Q is related to the 
cylinder self-energy, which will be important in obtain- 
ing a convergent energy expression for the system in the 
simulations (Section fV Bl and Appendix ID)| . 



B. Dimensionless description 

The parameter space of the system may be spanned by 
a minimal set of independent dimensionless parameters 
obtained from the ratios between characteristic length 
scales. These length scales are the rescaled Bjerrum 
length, q^in, the Gouy-Chapman length 
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(6) 



and the radius of the charged cylinder, R, and that of 
the outer boundary, D. The rescaled cylinder radius 



(7) 



equals the Manning parameter, ^. The ratio between the 
rescaled Bjerrum length and the Gouy-Chapman length, 
gives the so-called electrostatic coupling parameter 



(ij) 



2-Kq 



3«2 , 



'-B'^s, 



(8) 
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which can identify the importance of electrostatic corre- 
lations in a charged system [s^, IH, , and the ratio 
between D and i?, which enters only through the lateral 
extension parameter 



In 



(9) 



characterizing lateral finite-size effects. The relevant 
infinite-system-size limit is obtained for A ^ oo [T^ Isol 

EH. 

We shall use the dimensionless form of the Hamiltonian 
obtained by rescaling the spatial coordinates as x = x//i 
Isl. that is 



n 



N 



N 



s5]«3D(i,-x,) + 2e5]ln(4) +Co. (10) 



i=l 



The electroneutrality condition ^ in rescaled units reads 



2<L = 27rSiV, 



(11) 



where the left hand side is simply the rescaled area of 
the cylinder covered by the electric charge. The thermo- 
dynamic limit is obtained for N ~^ oo and L ^ oo, but 
keeping N/L = r/q (or equivalently, N/L — ^/S) fixed. 



C. CCT as a generic binding-unbinding process 

The statistical physical properties of the system may 
be investigated using the canonical partition function, 



Nl 



N 



: df,: f,; 



i=l 



exp 



n 



N 



(12) 



represented in cylindrical coordinates x^ — (fi, (jji, Zi), 
with the spatial integral running over the volume, V, of 
the space accessible for counterions, i.e. R < f < D. 

Naively, one may conjecture that the partition function 
(|12|l diverges in a certain range of Manning parameters, 
when the upper boundary of the radial integrals, D, tends 
to infinity, as may be indicated by the logarithmic form of 
the counterion-cylinder interaction, which gives rise to al- 
gebraic prefactors of the form r^ in the integrand. The 
possible emergence of a divergency in a charged cylindri- 
cal system was first pointed out by Onsager and the con- 
nection with the counterion-condensation transition was 
discussed by Manning 

Here we demonstrate this peculiar point using a trans- 
formation of coordinates, which provides the basis for our 
numerical simulations considered later in Sections IVl and 
IVII The radial coordinate is transformed as 



(13) 



upon which the partition function in 112|) transforms as 



3N td2N 



nT 



N 



nd 



exp 



N 



(14) 

where the volume integral runs over the region < y < 

A = \n{D/R), and 



n 



N 



N 

E 

1=1 



W{yi) -I- ^^t;3D(Xi - Xj 

iv) 



(15) 



is the transformed Hamiltonian of the system with 



(16) 



As seen, the original partition function is now mapped 
to the partition function of a system of interacting 
(repelling) particles in a linear potential well, W{y). 
This virtual potential includes the contributions asso- 
ciated with the cylindrical boundary, namely, the bare 
counterion-cylinder attraction (i.e. 2^y) and an entropic 
(repulsive) term from the measure of the radial integral 
(i.e. 2y), which may be regarded as an induced centrifu- 
gal component. 

For small Manning parameter, ^ < 1, the potential 
well, W{y), becomes purely repulsive suggesting that 
counterions unbind (or "de-condense") from the cen- 
tral cylinder departing to infinitely large distances as 
the outer confining boundary tends to infinity, A — 
\n{D/R) oo. In contrast for ^ > 1, the potential well 
exerts an attractive force upon counterions, which might 
lead to partial binding (or "condensation" ) of counterions 
even in the absence of confining walls. The new represen- 
tation of Zn in Eq. (|14|l . therefore, reflects the interplay 
between energetic and entropic factors on a microscopic 
level. 

Note that the rigorous analytical derivation of the 
aforementioned properties for counterions based on the 
full partition function is still an open problem, and only 
approximate limiting cases have been examined analyti- 
cally (Section HrEj. 



D. Onsager instability 

As a simple illustrative case, let us consider a "hypo- 
thetical" system, in which mutual countcrionic repulsions 
are switched off. The partition function H12|) thus factor- 
izes as Zn ^ Z^ , where 



Zi 



dye^'-'i^y= ' (17) 
z — 



is the single-particle partition function. It diverges for 
^ < 1, when the lateral extension parameter. A, tends to 
infinity, which implies complete de-condensation of coun- 
terions, i.e. the probability, P{r) ~ exp(— 2^1nr)/2^i, of 
finding counterions at any finite distance, r, from the 
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cylinder tends to zero (equivalent to a vanishing density 
profile, p{r) = NP{r)). But Zi and the counterionic 
density profile remain finite for ^ > 1, indicating that 
the Manning parameter = 1 is the onset of the CCT 
on the one-particle level, which we term here as the On- 
sager instability (in the spirit of Onsager's original argu- 
ment I'l ) . Onsager instability captures the basic features 
of the CCT. It exhibits the weak logarithmic convergence 
(via A = In I? /R) to the critical limit as the volume per 
polymer (~ D^) goes to infinity [G^I, and as shown in Ap- 
pendix ^ displays algebraic singularities in energy and 
heat capacity (at = 1) that may be identified by a set 
of scaling exponents. Such scaling relations are crucial in 
our analysis of the CCT in the following sections. 

We emphasize here that the results obtained within 
Onsager instability are by no means conclusive as soon 
as inter-counterionic interactions are switched on, which, 
as will be shown, lead to qualitative differences. In par- 
ticular, it turns out that a diverging partition function is 
not necessarily an indication of the onset of the CCT as 
asserted by the single-particle argument . 



cal phenomena |49j, can be shifted from its mean- field 
value due to inter-particle correlations for large cou- 
plings. Also it is interesting to examine whether the 
CCT exhibits scale-invariant properties near and if 
it can be classified in terms of a universal class of critical 
exponents. Such scaling relations are known to repre- 
sent relevant statistical characteristics of systems close 
to continuous phase transitions |49|. 

To address these issues, one has to define quantities 
which can serve as order parameters of the CCT. In the 
following section, we shall introduce such quantities and, 
by considering the mean-field theory, show that the or- 
der parameters indeed exhibit scaling behavior near the 
CCT critical point. We return to the influence of electro- 
static correlations on the critical Manning parameter and 
scaling exponents of the CCT in the subsequent sections. 

III. MEAN-FIELD THEORY FOR THE 
COUNTERION-CONDENSATION TRANSITION 

A. Non-linear Poisson-Boltzmann (PB) equation 



E. Beyond the Onsager instability: many-body 
effects and electrostatic correlations 



Many-body terms involved in the full partition func- 
tion (|12|l render the systematic analysis of the CCT 
quite difficult. The analytical results are available in the 
asymptotic limits of i) vanishingly small coupling param- 
eter 2^0, which leads to the mean-field or Poisson- 
Boltzmann (PB) theory, and ii) for infinitely large cou- 
pling parameter S — > cxd, which leads to the strong- 
coupling (SC) theory 52]. In the mean-field approxi- 
mation (case i), statistical correlations among countcri- 
ons are systematically neglected. In the opposite limit of 
strong couphng (case ii), the leading contribution to the 
partition function takes a very simple form comprising 
only the one-particle (counterion-cylinder) contributions, 
which is associated with strong electrostatic correlations 
(pronounced correlation hole) between counterions at the 
surface IH, |5^ |5^ [s^l • We study the mean-field pre- 
dictions for the CCT in Section ITTll The SC description 
(S oo) resembles the Onsager instability and will be 
discussed in Section Hvl and Appendix 1X1 The perturba- 
tive improvement of these two limiting theories in a sys- 
tem of finite coupling parameter, S, is formally possible 
by computing higher-order correction terms as previously 
performed for planar charged walls ^2)113) but will not 
be considered here. 

Interestingly, in both limits, the onset of the CCT is 
obtained as = 1, which is due to the simplified form 
of the counterionic correlations. An important ques- 
tion is whether the critical value, ^c, varies with the 
coupling parameter. Such a behavior may be expected 
since the Manning parameter represents the rescaled in- 
verse temperature of the system (i.e. ^ = T*/T with 

= qTe^ / {Aires ak-Q)) , which, as known from bulk criti- 



The mean-field theory can be derived systematically 
using a saddle-point analysis in the limit S ^ [s^- It 
is governed by the well-known Poisson-Boltzmann (PB) 
equation |50l l5l| . which, in rescaled units, reads (Ap- 
pendix 0) 

V| V 2CT(i) - fi(i) e-'^(*) (18) 
for the dimensionless potential field Here 

cr(i) = 5(f-^) (19) 

is the rescaled charge distribution of the cylinder and 

1 R<r < D, 



= n{r) 



(20) 



otherwise 



specifies the volume accessible to counterions. In the 
canonical ensemble, one has 



Y 



27rCL 



J dxr2(x) exp(— ■)/)) 



(21) 



Assuming the cylindrical symmetry (for an infinitely 
long cylinder) and using Eq. (|18|l and the global elec- 
troneutrality condition (|ll|l . one obtains 



f^) =2^ and (r'^^ 



0, 



(22) 



which are used to solve the PB equation (|18|l in the non- 
trivial region R<f<D^ i^. Th ereby, one obtains 
both the free energy (Section IIII C 2|l and the rescaled 
radial density profile of counterions around the charged 
cylinder 



-n{f)e 



(23) 
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The rescaled density profile, p(r), is related to the ac- 
tual number density of counterions, p{r), through p(f) — 
p(r)/(27r£Bcr2) ;52] (Appendixp. 

As shown by Alfrey et al. 50] and Fuoss et al. 'Fl'l , the 
PB solution takes different functional forms depending on 
whether ^ lies below or above the Alfrey- Fuoss threshold 



Therefore, for Manning parameter C < li one may use the 
first relation in Eq. H26(l to obtain the limiting behavior 
of the integration constant (3 as f Appendix IC l|l 



(31) 



A 



AF 



1 + A' 



(24) 



that is 



when A — > cx). Using this into Eq. 1)28(1 . one finds that 
the density of counterions at contact, ppb{R), asymptot- 
ically vanishes. Hence, the density profile (|23|l at any 
finite distance from the cylinder tends to zero for ^ < 1, 
i.e. 



ln[|^sinh2(;31n-|+coth-i^)] ^ < Aaf, 
V'PB(r) = ( (25) 
[ln[^sin2(/31ni-f cot-ii^)] ^ > Aaf, 

where f3 is given by the transcendental equations 

^/3coth(-/3A) ^ - A-AF, 

(26) 



PPBif 



0, 



(32) 



1+13^ 



l-/3cot(-/3A) 



AF- 



The PB density profile of counterions, Eq. ((23|l . is then 
obtained ior R < r < D as 



(3^ 



sinh-2(/31ni +coth-i ^) C < Aaf, 



(27) 



sin-2(/31ni+cot-ii^) C > Aaf, 



where we have arbitrarily chosen ippBi^ = i?) = to fix 
the reference of the potential. This condition also fixes 
K in Eq. I|25|l as well as the radial density of counterions 
at contact with the cylinder using Eq. (|23|l . i.e. 



_=ppB(i?) = -X 



(^_1)2_^2 ^<Aaf, 

(e -!)' + /?' e>AAF. 



(28) 



The density profiles given in Eq. (|27|l are in fact normal- 
ized to the total number of counterions, A^, a condition 
imposed via Eq. (|21|l . Using Eq. H23|) , the normalization 
condition in rescaled units reads f Appendix IB|| 



/ dffppB(r)=C- 
Jii 



(29) 



B. Onset of the CCT within mean-field theory 

The threshold of CCT within the mean-field PB theory 
was considered by several workers 9, 10, 11, 12, 13. 5CL 
Isij l . It may be obtained from the asymptotic behavior of 
the density profile (A oo) as reviewed below. 

First note that for A ^ 1, the Alfrey- Fuoss threshold 
Aaf, Eq- HH), tends to unity from below, i.e. 



A 



AF 



1 

A 



-0(A-2). 



(30) 



representing the de-condensation regime in the limit A — > 
oo. For C > 1, on the other hand, one has /3 — > for 
increasing A f Appendix IC l|l . and thus using Eq. ((28|l . 



ppB(i?) 



(33) 



Using the second relation in Eq. H27|) and expanding for 
small /3, the radial density profile follows as \wL l61| 



Ppb(^) 



(c-1)^ 




-2 r 


e 







l + (e-l)ln^ 
R 



(34) 



in the limit A oo (see also Appendix IC 4|) , which is 
finite and indicates condensation of counterions. This 
proves that the mean-field critical point is given by 



= 1, (35) 
corresponding to the mean-field critical temperature 



T, 



PB 



qre 



(36) 



C. Critical scaling-invariance: Mean-field 
exponents 

It is readily seen from Eqs. ^ and l|21l that 
the asymptotic density of counterions admits a scale- 
invariant or homogeneous form with respect to the re- 
duced Manning parameter, 



PB 



(37) 



close to the critical value — 1 . Note that the reduced 
Manning parameter equals the reduced temperature of the 
system, t = 1 — (T/T^^), when other quantities such as 
the dielectric constant, e, and the linear charge density of 
the cylinder, r, are kept fixed. (Experimentally, however, 
the Mannin g pa rameter may be varied by changing e [T^ 
|20| or T [l^ll8llT9ll20l | at constant temperature, in which 
case, ( can be related to the reduced dielectric constant 
or the reduced linear charge density.) 



7 



In a finite confining volume (finite A), such scaling 
forms with respect to ( do not hold since the true CCT 
is suppressed. Yet as a general trend we expect that 
for sufficiently large A, the reminiscence of such scaling 
relations appears in the form of finite-size-scaling rela- 
tions near the transition point. These relations would 
involve both C and the lateral extension parameter, A, 
(as the only relevant parameters in the mean-field limit) 
in a scale-invariant fashion as will be shown below. 



1. The CCT order parameters 

As possible candidates for the CCT "order parameter" , 
we use the inverse moments of the counterionic density 
profile 



1 



f df f " p(r) 
li rdrp(r) 



(38) 



where n > |^ . Note that these quantities reflect mean 
inverse localization length of counterions. In the con- 
densation phase (where counterions adopt a finite den- 
sity profile), one has Sn > 0, reflecting a finite local- 
ization length. But at the critical point and in the de- 
condensation phase (with vanishing counterionic density 
profile), one has S'„ = in the limit of infinite system 
size A — > oo, which indicates a diverging counterion lo- 
calization length. 

In order to derive the mean-field finite-size-scaling re- 
lations for Sn near = 1, we focus on the PB solution 
in the regime of Manning parameters ^ > Aaf, since for 
any finite A, we have Aaf < = 1 from Eq. (|Sn|l . 
Inserting the first relation in Eq. H27|) into Eq. (|38|l . we 
obtain 



«2 pD 

SP^ = L^I df f-"-^sm-^(j3\n^ +cot 

R 



£, Jr 



/3 



(39) 

Changing the integration variable as y = ln(f/i?), we get 



^PB 



n+l 



dye""^ sin~^(/3?/ cot 



-1^-1. 
/3 ' 



(40) 



For A ^ 1, the above relation may be approximated by 
a simple analytic expression as f Appendix IC 3() 



^r(c,A)-- 

n 



C' + /?'(C,A) 



(41) 



for ^ being sufhciently close to the critical value = 1. 

Using the above result, we may distinguish two limit- 
ing cases, where different scaling relations are obtained, 
namely, i) when A ^ oo but ^ = 1 — C^^/^ is finite and 
close to the critical value — 0, and ii) when A is fi- 
nite and large, but the system tends towards the critical 
point, C^Cc^ = 0. 



In the first stated before, we have /? — > for 

the above-threshold regime, C > 0; thus using Eq. H41|l . 
we obtain 



5r(C,A^(X3)c.^ 



(42) 



On the other hand, S^^^ vanishes for C 1^ (Appendix 
IC 3|l . Hence, the following scaling relation is obtained in 
the infinite-system-sizc limit A — > cxd. 



C^^/n 0<C<1, 







C<o, 



(43) 



which introduces the mean-field critical exponent asso- 
ciated with the reduced Manning parameter, C (or the 
reduced temperature, t) as 



XPB = 2. 



(44) 



The mean-field counterion-condensation transition is 
therefore characterized by a diverging (localization) 
length scale ^ C~^; as the critical point is ap- 

proached from above. The scaling relation H43II may also 
be derived in a direct way by considering a strictly infi- 
nite system (A = oo) as shown in Appendix IC 41 

In the limiting case (ii) with ( — > = 0, we have 
from Eq. ^ that (3 ~ 7r/(2A) when A is finite but 
large, A >• 1 f Appendix IC Therefore, using Eq. lUT)) 
we obtain 



5^(0, A) 



4nA2' 

which introduces a new scaling relation 

5P^(0,A) - A-''™ 
with the mean-field critical exponent 
7PB = 2 



(45) 



(46) 



(47) 



associated with the lateral extension parameter, A. This 
relation shows that the approach to the true CCT limit 
(when vanishes at the critical point) is logarithmi- 
cally weak as the box size, D, increases to infinity, i.e. 
5r(C = 0)^l/(ln i?/i?)2. 

The scaling relations (021 and indicate that S^^ 
takes a general scale-invariant form with respect to and 
A as 



^PB 



'VniCA 



7pb/xpe 



) 



(48) 



for sufficiently large A and in the vicinity of the mean- 
field critical point. The scaling function, I?„(it), has the 
following asymptotic behavior 



Vniu) 



const. M ^ 0, 



(49) 
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In general, the scale-invariant relations such as Eq. 
(I48|l may be obtained within the PB frame-work using 
the fact that the integration constant /3(C,A) takes a 
scale-invariant form as 



(50) 



Here B{u) is a scaling function which behaves asymptot- 
ically as (Appendix Q 



const, u 0, 



Bin) 



(51) 



Combining Eqs. H41|) and (|50|l . the scaling function 
T>n{u) is obtained in terms of B{u) as 



n '■ 



(52) 



which reproduces Eq. (|49|l when combined with Eq. (|51|l . 

The mean-field critical exponents xpb and 7pb appear 
to be independent of the order of the density moments, 
n. They may be used to characterize the mean-field uni- 
versality class of the CCT in all dimensions, since the PB 
results are independent of the space dimensionality. 



2. Mean-field energy and heat capacity 

As shown in a previous work , the mean- field canon- 
ical free energy of the counterion-cylinder system may 
be obtained using a saddle-point analysis from the field- 
theoretic representation of the partition function when 
S — !■ [521 ■ The rescaled PB free energy defined as 

netic energy part) 



IF^/lNk-oT), is given by (up to the trivial ki- 



^PB _ _ 1 
- In 



■ df 



PB 



df 



D 



f df e 



5{f - R^ippsif 



In 



(53) 



where Vcy\ = tiB?L is the actual volume of the cylinder. 
In the thermodynamic limit — > oo, the ratio Vcyi/N is 
a constant and will be dropped in what follows. 

Inserting the PB potential field, Eq. (|25|l . into the free 
energy expression (|53|l . we find that for ^ > Aaf 



-f-PB _ _ ^ 



1 



(1 -/32)A-f In 



+ ln[(e - 1)2 + /32] 
While for ^ < Aaf, we have 



1-1- /32 
ln(2e). 



(54) 



PB 



1 
1 

ln[(e 



(l+/32)A + ln 



(C-1) 



1 - 

ln(2e). 



(55) 



These expressions (up to some additive constants) have 
also been obtained by Lifson et al. js^ using a charging 
process method. 

The rescaled (internal) energy, E^^ = E^^/{NkBT), 
and the rescaled excess heat capacity, C^^ = 
/{NkB), can be calculated using the thermodynamic 
relations 



E 



PB 



PB 



^PB ^ _^ 



(56) 
(57) 



where the derivatives are taken at fixed volume, number 
of particles, and also for fixed charges and dielectric con- 
stant. A closed- form expression may be obtained for en- 
ergy using the relation E = (eeo/2) / dx (Vi/'cicc)^, where 
"^oicc = kBTippB/qe is the potential field in actual units. 
In rescaled units, the result is 



E 



1 

4? 



f df 



dv- 

df 



1 



(l-/32)A + ln 



(58) 

^ e<AAF, 

(59) 

e C>Aaf. 



In general, the above quantities can be calculated nu- 
merically using the transcendental equation H2t)|) . But in 
the limit of A — > oo, one may use the asymptotic results 
for (3 (Appendix [nj to derive the asymptotic form of the 
rescaled PB free energy as 



PB 



(^ - 2)A ^ < = 1, 
-A/e e > = 1- 



(60) 



The rescaled PB energy asymptotically behaves as 



^PB 



(61) 



and the rescaled PB excess heat capacity as 

e < ^r, 



(7™ = 



2A/^ e > i 



PB 



(62) 



The above results show that both energy and excess 
heat capacity develop a singular peak at the Manning 
parameter = i when the critical limit A — s- oo is ap- 
proached. The PB results also show that the free energy 
diverges with A both above and below the mean-field crit- 
ical point, in contrast with the behavior obtained within 
the (one-particle) Onsager instability , which suggests 
a connection between the onset of the counterion con- 
densation and the divergence of the partition function 
fSection Inland Appendix 0. 
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Another important point is that the PB heat capac- 
ity exhibits a discontinuity at = 1. Therefore, the 
CCT may be considered as a second-order transition as 
also pointed out in a previous mean- field study 38] . We 
shall return to the singular behavior of energy and heat 
capacity later in our numerical studies. 



IV. STRONG-COUPLING THEORY FOR THE 
CCT 

In the limit of large coupling parameter, S — > oo, 
the partition function of a charged system adopts an 
expansion in powers of the leading term of which 
comprises only single-particle contributions, i.e. a single 
counterion interacting with fixed charged objects |fi2Ll5.l| |. 
This leading-order theory, referred to as the asymptotic 
strong-coupling (SC) theory, describes the complemen- 
tary limit to the mean-field regime, S 3> 1, where inter- 
particle correlations are expected to become important 

The rescaled SC density profile for counterions is ob- 
tained as 1531 



(63) 



where u{r) = 2^ ln(f/i?) is the single-particle interaction 
energy and Aq is a normalization factor, which is fixed 
with the total number of counterions. Thus we have 



An = 



2(e-i) 



1 - e-2(€-i)A 



(64) 



in the cell model considered here. Note that for A — > cx3, 
Ao, and therefore the whole density profile, vanishes for 
^ < 1. But for ^ > 1, we get Aq 2(^ — 1)/^ and hence 
a finite limiting density profile as 



psc{f) 



2(g-l) 



-2« 



(65) 



This shows that the CCT is reproduced within the SC 
theory as well, and surprisingly, the critical value is found 
to be ^f*^ = 1 in coincidence with the mean-field predic- 
tion. Note however that the SC profile for ^ > 1 indicates 
a larger contact density for counterions as compared with 
the mean-field theory, e.g., for A ^ oo, one has 



PsciR) = 



(66) 



which is larger than the PB value (|33|l by a factor of 

Psc{R)/ppBiR) = 2(1 - The SC density profile 

also decays faster than the PB profile indicating a more 
compact counterionic layer at the cylinder. This reflects 
strong ionic correlations in the condensation phase (f > 
1) for S ^ 1 as will be discussed further in the numerical 
studies below. 



Using Eq. (|63|l . the SC order parameters can be cal- 
culated as 

2(^ _ 1) 1 - g-(2C-2+n)A 

(e, A) = _2 + n) 1 - e-2(«-i)A (^^^ 

for arbitrary ^ and A. For A — s- oo, vanishes for 



^^''(Coo) 



2(^-1) 



e(2e 



(68) 



for ^ > (^f^ = 1. In the vicinity of the critical point, S^'^ 
exhibits the scaling relation 



5^(C,oo) 



2C 



(69) 



which gives the SC critical exponent associated with the 
reduced Manning parameter, C. — I — Cc/C; Xsc — 1- 
In a finite system and right at the critical point ~ ^f*-", 
S^'^ exhibits the finite-size-scaling relation 



^^(0,A).-, 



(70) 



which gives the SC critical exponent associated with the 
lateral extension parameter. A, as 7sc = 1. 

These exponents are different from the corresponding 
mean-field values, Eqs. (|44|l and H47|l . and in fact coincide 
with the Onsager instability results. As will be shown 
later, the SC predictions in fact break down near the 
CCT critical point. 



V. MONTE-CARLO STUDY OF THE CCT IN 

3D 

The preceding analysis within the mean-field and the 
strong-coupling theory reveals a set of new scaling rela- 
tions associated with the counterion-condensation tran- 
sition (CCT) in the limit of infinitely large (lateral) sys- 
tem size. In the following sections, we shall use numer- 
ical methods to examine the critical behavior in various 
regimes of the coupling parameter, and thereby, to exam- 
ine the validity of the aforementioned analytical results. 



A. The centrifugal sampling method 

The major difficulty in studying the CCT numerically 
goes back to the lack of an efficient sampling technique. 
Poor sampling problem arises for counterions at charged 
curved surfaces in the infinite-confinement-volume limit 
because, contrary to charged plates, a finite fraction of 
counterions always tends to unbind from curved bound- 
aries and diffuse to infinity as the system relaxes to- 
ward its equilibrium state. This situation is, of course, 
not tractable in numerical simulations; hence to achieve 
proper equilibration within a reasonable time, charged 
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cylinders are customarily considered in a confining box 
(in lateral directions) of p ractically large size. As well 
known H [H El 11 113, lateral finite-size effects 
are quite small for sufficiently large Manning parameter 
> ^c)- But at small Manning parameters ~ ^c), 
these effects become significant and suppress the de- 
condensation of counterions. 

The mean-field results already reveal a very weak 
asymptotic convergence to the critical transition con- 
trolled by the logarithmic size of the confining box A — 
\n{D/R). Hence one needs to consider a confinement vol- 
ume of extremely large radius, D, to establish the large-A 
regime, where the scaling (and possibly universal) prop- 
erties of the CCT emerge. For this purpose, clearly, the 
simple-sampling methods within Monte-Carlo or Molec- 
ular Dynamics schemes l23i, 121 HI HI [gl |65, £6J are 
not useful at low Manning parameter as they render an 
infinitely long relaxation time. 

We shall therefore introduce a novel sampling method 
within the Monte-Carlo scheme, which enables one to 
properly span the relevant parts of the phase space for 
large confining volumes. In three dimensions, we use the 
configurational Hamiltonian l|10|) in rescaled coordinates. 
The sampling method, which we refer to as the centrifu- 
gal sampling^ is obtained by mapping the radial coordi- 
nate to a logarithmic scale according to Eq. (|13|l . i.e. 
y — ln(f/i?), which leads to the transformed partition 
function H14I) . As explained before (Section III C|l . the 
entropic (centrifugal) factor, exp(2y), is absorbed from 
the measure of the radial integrals into the Hamiltonian, 
yielding the transformed Hamiltonian in Eq. (|15|l . 

We thus simulate the system using Metropolis algo- 
rithm 67] , but making use of the transformed Hamilto- 
nian IjlSfl . The entropic factors, which cause unbinding 
of counterions, are hence incorporated into the transition 
probabilities of the associated Markov chain of states, 
that generates equilibrium states with the distribution 
function ~ exp(— /37i^). The averaged quantities, say 
Q, follow by extracting a set of T values {Qi, . . . , Qt} in 
the course of the simulations as Q = X^tLi Qt/T, which, 
for sufficiently large T, produces the desired ensemble 
average (Q), i.e. 



Q 



1 ^ 



3N rf2N 



N 



N\Z_ 



N 



(Q>, 



(71) 



up to relative corrections of the order l/\/T 



B. Simulation model and parameters 

The geometry of the counterion-cylinder system in our 
simulations is similar to what we have sketched in Figure 
n We use typically between A'' = 25 to 300 counterions 
(most of the results are obtained with A^=100 and 200 
particles) and increase the lateral extension parameter. 



A = ln(D/i?), up to A = 300. We also vary the Manning 
parameter, ^, and consider a wide range of values for the 
electrostatic coupling parameter, 2, from S = 0.1 (close 
to the mean-field regime) up to S = 10^ (close to the 
strong-coupling regime). 

The cylindrical simulation box has a finite height, L, 
which is set by the electroneutrality condition Hll|l . i.e. 
L — N'E./^. In order to mimic the thermodynamic limit 
and reduce the finite-size effects due to the finiteness of 
the cylinder height, we apply periodic boundary condi- 
tions in z direction (parallel to the cylinder axis) by repli- 
cating the main simulation box infinitely many times in 
that direction. The long-range character of the Coulomb 
interaction in such a periodic system leads to summation 
of infinite series over all periodic images. These series are 
not generally convergent, but in an electroneutral system, 
the divergencies cancel and the series can be converted to 
fast- converging series. We use the summation techniques 
due to Lekner [b^I and Sperb [i^, which are utihzed to 
the one-dimensionally periodic system considered here- 
see Appendix ^ (similar methods are developed in Ref. 
[70|). Finally in order to obtain reliable values for the 
error-bars, the standard block-averaging scheme is used 
[tJ. The simulations typically run for ~ 1.1 x 10^ Monte- 
Carlo steps per particle with ~ 10^ steps used for the 
equilibration purpose. 



VI. SIMULATION RESULTS IN 3D 
A. Overall behavior in the inflnite-system-size limit 

1. Distribution of counterions 

Let us start with the distribution of counterions as 
generated by the centrifugal sampling method. In Fig- 
ure 121 typical simulation snapshots are shown together 
with the countcrionic density profile for small coupling 
parameter S = 0.1. Counterion distribution is shown 
for large (A = 100), intermediate (A = 25) and small 
(A — 10) lateral extension parameter. The counterion- 
condensation transition is clearly reproduced for large A 
(Figure [2^-): counterions are "de-condensed" and gather 
at the outer boundary at small Manning parameter 
(shown for ^ — 0.7), while they partially "condense" and 
accumulate near the cylinder surface for large Manning 
parameter (shown for ^ = 2). The Manning parameter 
^ = 1, as seen, represents an intermediate situation. This 
trend is demonstrated on a quantitative level by the ra- 
dial density profile of counterions p{f) (Figure [5J:, main 
set), which tends to zero by decreasing ^ down to about 
unity. Note that relatively large fiuctuations occur at low 
^ making p{f) an inconvenient quantity to precisely lo- 
cate the critical value ^c, which will be considered later. 
The data moreover follow the mean-field PB density pre- 
diction, Eq. lf77|l . shown by solid curves, as expected 
since the chosen coupling parameter is small. 
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r/R 

FIG. 2: Typical snapshots from the simulations on the counterion-cylinder system in 3D for a) lateral extension parameter 
A = \n{D / R) = 100 and three different Manning parameters ^ = 0.7, 1.0 and 2.0 as indicated on the graph, and b) for Manning 
parameter ^ = 1.0 and smaller lateral extension parameters A = 10 and 25. The snapshots show top-views of the simulation 
box (see Section IV Bl and Figure with radial distances shown in logarithmic scale y — ln{f/R). Point- like counterions are 
shown by black spheres and the charged cylinder by a circle in the middle. Figure c) gives the simulated radial density of 
counterions in rescaled units, p(f) = p(r)/(27ri'BO"g ), as a function of the (linear) distance from the cylinder axis. Main set 
shows the data for A = 100 and Manning parameters ^ — 2.0 (open square), 1.5 (open triangle-ups) , 1.1 (open diamonds), 1.0 
(filled squares) and 0.7 (filled circles) from top to bottom. Inset shows the same for ^ = 1.0, but for various lateral extension 
parameters A =10, 25 and 100 (top to bottom). Solid curves represent the mean-field PB prediction Eq. H27|l . Number of 
counterions here is N = 100 and the coupling parameter H = 0.1. Error-bars are smaller than the size of symbols. 



The transition regime at intermediate ^ exhibits strong 
finite-size effects. As may be seen from the snapshots 
in Figure HJd, the de-condensation process at ^ = 1 
is strongly suppressed for smaU logarithmic sizes A = 
\n{D/R) — 10 and 25. The corresponding density pro- 
files (inset of Figure ^) indicate a sizable accumulation 
of counterions near the cylinder surface, which is washed 
away only by taking a sufficiently large A. Such finite- 
size effects at low ^ are also observed in previous numeri- 
cal studies, which have devised simulations in linear scale 
and thus considered only small confinement volumes per 
polymer (typically A < 10) [H IH HI HI HI . In some 
studies j72ij , these effects have been interpreted as an ev- 
idence for counterion condensation at small ^, leading to 
the incorrect conclusion that no condensation transition 
exists. 



2. Condensed fraction of counterions 

Our results for large ^ exhibit a counterionic density 
profile that extends continuously from the cylinder sur- 



face to larger distances. This indicates that making 
a distinction between two layers of condensed and de- 
condensed counterions, in the sense of two-fiuid models 
freouently used in literature HIllllllEIEIillllllll 
mllliiEilillill, requires a criterion. 

The two-fluid description predicts a fraction of 

fO ^ < 1 

au = < (72) 

ii-i/e e>i 

of counterions to reside in the condensed layer (which is 
considered as a layer with small thickness at the poly- 
mer surface), when the infinite-dilution limit is reached. 
Previous studies .9, IQu JJ,^ 13, 24, 43] show that the 
Manning condensed fraction, am, may also be identified 
systematically within the Poisson-Boltzmann theory by 
employing an inflection-point criterion l4^ . This can 
be demonstrated using the PB cumulative density (the 
number of counterions inside a cylindrical region of ra- 
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FIG. 3: Cumulative density, n(f), per total number TV, of 
counterions as a function of the logarithmic distance from the 
charged cylinder, ln(f/i?). Dot-dashed curves are simulation 
results for H = 0.1, TV = 70, and A = 300 and for various 
Manning parameters as shown on the graph. These curves 
also closely represent the PB prediction I73II . which are not 
explicitly shown. Solid curves show the simulation data for 
large coupling parameter H = 10^ and for ^ = 3.0 and ^ = 2.0. 
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FIG. 4: Main set: Simulated condensed fraction of coun- 
terions, Q, as defined via the inflection-point criterion, as a 
function of Manning parameter, ^, for H = 0.1 (diamonds). 
Solid curve displays Manning (or the PB) limiting value, qm , 
for A = ln(D/i?) -> oo (Eq. IJU). Inset: Condensed fraction 
as a function of the coupling parameter, H, for ^ = 3.0. These 
data are obtained for A = 300 and TV = 70. 



dius r), npB(r), obtained as 
npB(v) 27rL 



N 



r' dr' ppB{r') 



(C- 1) -/3coth 
(C- 1) -/3cot 



l3y + coth ^ ^ 



I3y + cot-' ^ 



e < Aaf, 
(73) 
e > Aaf, 



using Eq. (|27|l . For ^ > Aaf, npB{r) exhibits an inflec- 
tion point at a radial distance when plotted as a func- 
tion of y = ln(r/i?) (2^. One can show that for A oo, 
only the fraction of counterions, that lie within the cylin- 
drical region r < , remains associated with the cylinder 
and tends to the Manning condensed fraction, i.e. 



N 



(74) 



In other words, only this fraction of counterions con- 
tribute to the asymptotic density profile and the rest (1/^ 
of all) is pushed to infinity (Appendix IC 4|l . 

The simulations results for the cumulative density as 
a function of the logarithmic radial distance y = ln{f/R) 
are shown in Figure |31 for various Manning parameters 
(solid and dot-dashed curves). Here we have chosen a 
very large lateral extension parameter A = ln(_D/i?) — 
300, which exhibits the concept of condensed fraction 
more clearly. The data show an inflection point, which is 
located approximately at j/* = ln(r,/i?) ~ A/2 for large 
^ (for small ^ — » 1, the location of the inflection point, 

, tends to R-see Appendix IC 2\i . The rapid increase 
of n{f) at small (r ~ R) and at large distances (r ~ 
D) reflects the two counterion-populated regions at the 



inner and outer boundaries, which are separated by an 
extended plateau (compare with Figure For small A, 
the inflection point has a non- vanishing slope and the two 
regions are not quite separated (data not shown) [3, Q0| 
(see Ref. |35( for a similar trend in an extended two-fluid 
model) . 

Using the inflection-point criterion, the condensed frac- 
tion, a, may be estimated as a = n(r^)/N "2?!, which 
roughly corresponds to the plateau level in Figurc|21 Sim- 
ulation results are shown in Figure 01 for large A 300. 
Let us first consider the case of a small coupling param- 
eter S = 0.1, where the simulated cumulative density, 
n{f) (dot-dashed curves in Figure closely follows the 
PB prediction H73|l (PB curves are not explicitly shown). 
The calculated condensed fraction (diamonds in Figure 
0|| agrees already quite well (within < 1%) with the Man- 
ning or PB limiting value om (solid curve in Figure 01). 

An important question is whether the form of the cu- 
mulative density profile, n{f), and the condensed frac- 
tion are influenced by electrostatic correlations for in- 
creasing coupling parameter S. In Figure |21 we show 
n{f) from the simulations for S = 10^ and for two values 
of Manning parameter ^ = 2.0 and 3.0 (solid curves). 
This coupling strength generally falls into the strong- 
coupling regime for charged systems, where electrostatic 
correlations are expected to matter |5g (note that DNA 
with trivalent counterions represents S ~ 10'^, but with 
a larger ^ ^ 12). As seen, n{r) shows a more rapid 
increase at small distances from the cylinder (condensed 
region) indicating a stronger accumulation of counterions 
near the surface. This trend is also observed in previous 
simulations p3 . l64l l65|.|66ll and in experiments with mul- 
tivalent counterions |73j |. and will be analyzed in more 
detail in the following sections. 
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FIG. 5: Simulation data for the order parameter Si = {1/r) 
as a function of Manning parameter, ^, for various coupling 
parameters S — 0.1 up to 10^ as indicated on the graph. 
The mean-field PB and the strong-coupling predictions are 
calculated from Eqs. 1401 and 16 711 (solid and dashed curves 
respectively) . The lateral extension parameter is A = 300 and 
the number of counterions is A*' = 200 for H — 0.1, A*' = 50 
for H = 10^, and N = 70 for other coupling parameters. 
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FIG. 6: Radial density of counterions in rescaled units, 
p(f) — p(r)/(27r^BO"g ), as a function of the (linear) distance 
from the cylinder axis for Manning parameter ^ = 3.0 and 
various coupling parameters (E = 0.1 up to 10^) as shown 
on the graph. Here A — 300 and the number of counteri- 
ons is A'' = 50 for H — 10^, and A*' = 70 for other values 
of H. The mean-field (solid curve) and the strong-coupling 
(dashed curve) predictions are obtained from Eqs. 1271 and 
15311 ■ which, for A — 300, roughly coincide with the asymp- 
totic expressions 1341 and 1651 . 



However, in contrast to previous conclusions (obtained 
based on small values of A) 0jE9j the aforementioned 
behavior for large S does not imply a larger condensed 
fraction as defined within the inflection-point criterion. 
Since as seen in Figure |21 the large- distance behavior of 
the density profile is not influenced by electrostatic corre- 
lations, and so remains the condensed fraction (plateau 
level) unaffected for increasing coupling strength (inset 
of Figure 0}. This result can be appreciated only when 
the asymptotic behavior for A ^ 1 is considered. 

3. The order parameters 5„ 

The n-th-order inverse moment of the counterionic 
density profile may be calculated numerically using 

^« = ]^E-T" (75) 

4=1 

for n > 0, where fi is the radial distance of the z-th coun- 
terion from the cylinder axis and the bar sign denotes 
the Monte-Carlo time average after proper equilibration 
of the system. The overall behavior is shown in Figure 
13 for Si as a function of Manning parameter, ^. Re- 
call that a vanishing order parameter. Si, indicates the 
complete de-condensation of counterions, while a finite 
Si reflects a finite degree of counterion binding to the 
charged cylinder (corresponding to a finite localization 
length ~ 1/Si). 

As seen from the figure, de-condensation can occur in 
all relevant regimes of the coupling parameter S. For 
large Manning parameter, electrostatic coupling effects 



become important and shift the order parameter to larger 
values exhibiting a crossover from the mean-field pre- 
diction (solid curve), which is verified for small S < 1, 
to the strong-couplin g pr ediction (dashed curve) at very 
large values of S [s^ l53l l55l | . The mean- field result fol- 
lows from Eq. H4()|l and the strong-coupling prediction 
is obtained using Eq. (|67|) . As seen, in the transition 
regime ^ '--^ 1, the order parameter data remain close 
to the mean-field curve and deviate from the SC predic- 
tion. A close examination of correlation effects as well 
as finite-size effects in this region is quite important in 
determining the scaling behavior and will be considered 
later. Here we concentrate on the correlation-induced 
crossover behavior in the condensation phase. 



4- Electrostatic correlations at surface and for large ^ 

In Figure ini we plot the simulated radial density pro- 
file of counterions, p(f), for ^ — 3.0 and consider several 
different coupling parameters. In agreement with the pre- 
ceding results, the counterionic density in the immediate 
vicinity of the charged cylinder increases for increasing S 
exhibiting large deviations from the mean-field prediction 
(seeRef. [zl for a similar trend at charged plates). For 
a given surface charge density cts, the observed trend is 
predicted, e.g., for increasing counterion valency, q, since 
the coupling parameter scales as S ~ (Eq. JSJ). The 
crossover from the mean-field PB prediction (solid curve) 
to the strong-coupling prediction (dashed curve) appears 
to be quite weak, in agreement with the situation ob- 
served for counterions at planar charged walls [s^ . These 
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limiting profiles are calculated from Eqs. (|27|l and H63() 
respectively, and both exhibit an algebraic decay with the 
radial distance, f. But the SC profile shows a faster de- 
cay and thus a more compact counterion layer near the 
surface at large coupling strength (compare Eqs. H34() 
and jnSIl)- 

An interesting point is that the simulated density at 
contact with the cylinder shows a more rapid increase 
when the coupling parameter increases from S = 10 to 
S = 100 as compared with other ranges of S (Figure IH)). 
This is in fact accompanied by the formation of corre- 
lation holes around counterions near the surface as we 
show now. 

In order to examine countcrion-counterion correlations 
at the surface, we consider the one-dimensional pair 
distribution of counterions, giD(z), which measures the 
probability of finding two counterions lined-up along z- 
axis (i.e. along the cylinder axis with equal azimuthal 
angles (p) at a distance z from each other. In Figure [7| 
we plot the unnormalized pair distribution function de- 
fined via 

giD{S) = -^Yl y^^' ~ ^3 - ^) -<Pj)\ (76) 

where the prime mark indicates that the sum runs only 
over counterions at the surface (defined in the simula- 
tions as counterions residing in a shell of thickness about 
the Gouy-Chapman length, /i, around the cylinder). At 
small coupling parameter (S — 10, cross symbols), the 
pair distribution function only shows a very weak deple- 
tion zone at small distances along the cylinder axis. For 
larger values of S, one observes a pronounced correlation 
hole at small distances around counterions, where the 
distribution function vanishes over a finite range. This 
correlation hole develops in the range of coupling param- 
eters 10 < S < 100, which marks the crossover regime 
between the mean-field and the strong-coupli ng r egime 
(compare cross symbols and filled triangle-ups) |53| ■ The 
correlation hole appears only for sufficiently large Man- 
ning parameter ^ (large enough number of condensed 
counterions) and is distinguishable in our simulations for 
e > 1.2. 

The small-separation correlation hole is followed by 
an oscillatory behavior for elevated ^ indicative of a 
short-ranged liquid-like order among counterions line-up 
along the cylinder axis (distinguishable from the data for 

> 2.0 in the large-couphng regime S > 100). The loca- 
tion of the first peak of giD gives a rough measure of the 
typical distance between lined-up counterions, Oz, at the 
cylinder surface. This distance is set by the local elec- 
troneutrality condition a^r = q. In rescaled units, we 
obtain 



from Eqs. (|3J|, l|njl and ijHJl, which is used to rescale the 
horizontal axis of the graph in Figure [7| 



0.01 




FIG. 7: The one-dimensional pair distribution function of 
counterions at contact with cylinder as defined in Eq. 1)76^ . 
Symbols show simulation data for H = 10^ and ^ — 4.0 (filled 
diamonds), H = 10^ and ^ = 2.0,3.0 and 4.0 (open symbols 
from bottom to top), and for coupling parameters H = 10 and 
H — 100 with Manning parameter ^ — 3.0 (cross symbols and 
filled triangle-ups respectively). 



Note that the correlation hole size increases with the 
coupling parameter and thus counterions at the sur- 
face become highly isolated, reflecting dominate single- 
particle contributions for S ^ 1 ISZlSSl. In fact, as dis- 
cussed elsewhere jS^, EE , the single-particle 
form of the SC theory (obtained formally for S — > oo) is 
a direct consequence of large correlation hole size around 
counterions at the surface. Clearly, for the counterion- 
cylinder system, this can be the case only for sufficiently 
large Manning parameter, where a sizable fraction of 
counterions can gather near the surface. Consequently 
in this regime, the data tend to the strong-coupling pre- 
dictions for elevated S (Figures|Sland|ni as also verified in 
the simulations of charged plates, where all counterions 
are bound to the surface |^ , and two charged cylinders 
with large ^ [s^. This also explains why the SC theory, 
though being able to reproduce the CCT on a qualitative 
level, fails to capture the quantitative features near the 
critical point (except for the value of the critical Manning 
parameter), where counterion are mostly de-condensed. 

B. Critical Manning parameter 

We now turn our attention to the behavior of counte- 
rions near the critical point and begin with determining 
the precise location of the critical Manning parameter. 

To this end, we shall employ a procedure similar to the 
method of locating the transition temperature in bulk 
critical phenomena [t^- Namely, one expects that 
the transition point is reflected by a singular behavior 
in thermodynamic quantities such as energy or heat ca- 
pacity as already indicated by the mean-field results ob- 
tained in Section llTrn2l The mean (internal) energy, E, 
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FIG. 8: a) The rescaled internal energy, E — En /(NkBT) and b) the rescaled excess heat capacity, C — Cn /{NkB), of the 
counterion-cyhnder system as a function of Manning parameter, ^. Open symbols show the data for small coupling parameter 
S = 0.1 and for increasing lateral extension parameters A = 100,200 and 300 as indicated on the graph. Filled symbols are 
the data for large coupling H = 10^ and A = 200. Here number of counterions is A'' = 100. Solid curves show the asymptotic 
PB prediction for A — > oo, Eqs. I61II and 1621 . Dashed curves are the full PB result for A = 100 and 300 (from bottom to 
top), which are calculated numerically using Eqs. I59II and 18111 . The inset shows a closer view of the energy peak. 



and the excess heat capacity, C, may be obtained directly 
from the simulations and in rescaled units as 



namic relation 



E 



C = 



E 



N 



Nk^T 



n 



N 



NkB 



= N 



NksT 



\NkBT 



(78) 
(79) 



where the configurational Hamiltonian Tijv is defined 
through Eq. JTHl and STIn ='Hn - {Un)- 

Simulation results for the rescaled energy, E, and the 
rescaled excess heat capacity, C, in Figure |H1 (symbols) 
show a non-monotonic behavior as a function of ^. The 
energy develops a pronounced peak and the heat capac- 
ity exhibits a jump at intermediate Manning parameters, 
which become singular for A increasing to infinity. The 
general behavior of energy and heat capacity can be un- 
derstood using simple arguments as follows. 

For sufficiently small ^, counterions are all unbound 
and the electrostatic potential in space is roughly given 
by the bare potential of the charged cylinder, i.e. V'(^) — 
2^1n(f/i?). This yields the rescaled internal energy, E, 
(via integrating over the square electric field, Eq. 
as 



^.1 
4^ 



D 



f df ( — ^ 
dr 



(80) 



for A — \n{D/R) ^ 1. Intuitively, this result may be ob- 
tained also by assuming that counterions experience the 
potential of the cylinder at the outer boundary; thus one 
simply has E ~ 'tp{D)/2 ~ ^A, which explains the linear 
increase of the left tail of the energy curve with both ^ 
and A (Figure |S^). Now using the following thermody- 



dE 



E-C, 



(81) 



the excess heat capacity is obtained to vanish in the de- 
condensation regime, i.e. C* ~ (Figure IHti). Hence, the 
heat capacity reduces to that of an ideal gas of particles. 

For large ^, the electrostatic potential of the cylinder is 
screened due to counterionic binding. If we estimate the 
screened electrostatic potential of the cylinder as i^i^) — 
2 ln(f/i?), which can be verified systematically within the 
PB theory |^ , we obtain the energy and the heat 
capacity as 



E~A/^ and C ~ 2A/f . 



(82) 



These results may also be obtained by noting that only 
the fraction 1/^ of de-condensed counterions (Section 
IVI A 2|l contributes to the energy on the leading order; 
thus E ~ ^{D)/{2Q ~ A/^. The above asymptotic esti- 
mates in fact coincide with the asymptotic (A — > oo) PB 
results (|61|l and H62() , which are shown by solid curves in 
Figure |H1 

The preceding considerations demonstrate that the 
non-monotonic behavior of the energy and excess heat 
capacity results directly from the screening effect due to 
the condensation of counterions for increasing ^. Hence, 
the singular peaks emerging in both quantities reflect 
the onset of the counterion-condensation transition, ^c, 
which occurs in the thermodynamic infinite-system-size 
limit N ^ oo and A ^ oo. Within the PB theory (solid 
and dashed curves in Figure the location of the peak 
of energy, ^f''^^(iV, A), tend to the mean-field critical 
value = 1 from below as A increases obeying the 
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FIG. 9: Main set: Location of the peak of energy, S^f , as a 
function of the inverse lateral extension parameter, 1/A = 
l/ln(_D/_R). Open symbols are simulation results for small 
coupling parameter H = 0.1 and for various number of counte- 
rions (from bottom): A'^ = 25 (triangle-downs), 50 (squares), 
100 (triangle- ups), 200 (circles) and 300 (diamonds). Solid 
curve shows the mean-field PB prediction for the peak loca- 
tion as obtained by numerical evaluation of the full PB energy, 
Eq. 1591 . Inset: Location of the energy peak as a function of 
the coupling parameter, H, for = 200 and A = 300. 



finitc-size-scalmg relation 



PB 



(83) 



which is obtained using the full PB energy H59|l . On 
the other hand, the location of the peak of the PB heat 
capacity approaches from above. 

We locate the critical point from the asymptotic be- 
havior of the energy peak, , for increasing N and A. 
(The heat capacity peak is found to be located further 
away from the critical point than the energy peak, resem- 
bling the well-known behavior of the heat capacity peak 
in finite Magnetic systems [t^, which makes it inconve- 
nient for our purpose). In Figure |H1 we show the simu- 
lation results for (symbols) as a function of A^^ for 
S = 0.1 and for various number of particles (main set). 
These data are obtained using the thermodynamic rela- 
tion H81() , which allows us to calculate the first derivative 
of energy, dE/dS^, directly from the energy and the heat 
capacity data without referring to numerical differentia- 
tion methods which typically generate large errors near 
the peak. As seen, for increasing iV, the data converge to 
and closely follow the mean-field prediction (solid curve) 
within the estimated error-bars; for N > 100, lies 
within about 1% of the PB critical Manning parame- 
ter = 1. Since in the simulations we have used 
A < 300, the behavior of S^f for very small A~^ 
is not obtained, nevertheless, the excellent convergence 
of the data for S = 0.1 to the PB prediction gives an 
accurate estimate for the critical Manning parameter as 



= 1.00 ±0.002. 



(84) 



Our results for larger values of the coupling parame- 
ter, 2, in the inset of Figure El show that the location 
of the energy peak does not vary with the coupling pa- 
rameter. Therefore, we find that the critical Manning 
parameter is universal and given by the mean-field value 
= 1-0. Recall that the same threshold is obtained 
within the Onsager instability and the strong-coupling 
analysis (Sections III El and lIVp . 

Another important result is that the OCT is not as- 
sociated with a diverging singularity, in contrast to the 
Onsager instability prediction 1]. But, the energy at any 
finite value of ^, and also the heat capacity for ^ > 1, tend 
to infinity (as A) when the lateral extension parame- 
ter, A, increases to infinity, which, as illustrated before, 
refiects the logarithmic divergency of the effective elec- 
trostatic potential in a charged cylindrical system. The 
CCT, however, exhibits a universal discontinuous jump 
for the excess heat capacity at ^c, and thus indicates a 
second-order phase transition (Figure |HJ). 



C. Scale-invariance near the critical point 

Now that the precise location of the critical Manning 
parameter is determined, a finite-size analysis, similar to 
what we presented within the mean-field theory, may be 
used to determine the near-threshold properties of the 
CCT order parameters from the simulation data. 

Note that in the simulations, finite size effects arise 
both from the finiteness of the system size (via the lat- 
eral extension parameter. A), and also from the finite- 
ness of the number of counterions, N] the latter being 
related to the finiteness of the height of the main sim- 
ulation box L = Nq/r (Sectional, which has a sizable 
influence on the transition, although the implemented 
periodic boundary condition already reduces its effects. 
In what follows, we present the numerical evidence for 
scaling relations with respect to both N and A. The 
asymptotic behavior for increasing N and A to infinity 
provides us with the scaling behavior with respect to the 
reduced Manning parameter, C (or the reduced temper- 
ature, t), which characterizes the CCT universality class 
in 3D. 



1. Finite-size effects near 

In Figure [TUI (main set), we show the order parame- 
ter 5*1 as a function of 1/A and in the vicinity of the 
critical point = 1 (number of counterions N = 100 is 
fixed). 5*1, which represents the mean inverse localiza- 
tion length of counterions, gradually decreases with de- 
creasing 1/A as de-condensation of counterions becomes 
gradually more pronounced, but for Manning parame- 
ters as large as ^ = 1.05 (open circles), the data quickly 
saturate to a finite value as A oo. For sufficiently 
small Manning parameter (e.g., ^ < 0.97), on the other 
hand, 5*1 converges to zero. In the vicinity of the critical 
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FIG. 10: Main set: Order parameter Si — (l/r) as a function 
of the inverse lateral extension parameter 1/A. Open symbols 
are data for H — 0.1 and Manning parameters (from top): 
^ — 1.05 (squares), 1.01 (circles), 1.0 (open diamonds), 0.99 
(triangle- ups) and 0.97 (triangle-downs). Filled diamonds are 
the data for large coupling parameter H = 10^ and ^ = 1.0. 
Number of counterions = 100 is fixed. Dot-dashed curves 
are the PB result, Eq. 1401 . for the corresponding f . Inset: 5*1 
as a function of the inverse number of counterions 1/A^ for ^ = 
1.0 and for the coupling parameters H = 0.1 (open diamonds), 
10^ (filled diamonds) and 10^ (open circles). Dashed line 
shows the power-law exponent 2/3. Here A — 300. 



point = 1, diamonds), a non-saturating behavior is 
found suggesting a power-law decay as Si ~ A~^, where 
7 > 0. As seen, the data at ^ = 1 roughly coincide 
for both small coupling (S — 0.1, open diamonds) and 
large coupling (S = 10^, filled diamonds) indicating that 
electrostatic correlations do not influence the scaling be- 
havior (see below). There still remain non- negligible de- 
viations between the simulation data at the critical point 
(diamonds) and the PB power-law prediction H45I) with 
7pb — 2, which is shown in the figure by a straight dot- 
dashed line. These deviations arise from the finitcncss of 
the number of particles. 



Interestingly, the data obtained for increasing num- 
ber of counterions, N (at fixed lateral extension parame- 
ter, A), also indicate a power-law decay near the critical 
point, i.e. as 5*1 ~ , where v > Q. This is shown in 
the inset of Figure [TUI where the scaling exponent v ap- 
pears to be about 2/3 (represented by a dashed line). In 
fact, for sufficiently large N , the data deviate from this 
power-law behavior since finite-size effects due to lateral 
extension of the system. A, are simultaneously present. 
Thus in order to determine the exponents 7 and v, a 
more systematic approach is required, which should in- 
corporate both lateral-size and ion- number effects. 



2. Generalized finite-size scaling relations 

In brief, our data suggest that at the critical point 
(C = 1 — ^c/^ — 0) and for a bounded system (finite A) 
in the thermodynamic limit N 00, the order parameter 
5„(C,A,iV) = (1/f") decays as 



5„(0,A,oo)~ A- 

while in an unbounded system (A ^ 
N, we expect a power-law decay as 

5„(0,oo,iV) - N' 



(85) 

and for finite 



(86) 



In thermodynamic infinite-system-size limit (A — > 00 and 
N 00), the true critical transition sets in with 5„(C < 
0, 00, 00) = 0, and we anticipate the scaling behavior with 
the reduced Manning ( as 



£•„((, 00, 00) ^ C 



(87) 



in a sufficiently small neighborhood above — I- 

These scaling relations may all be deduced from a 
general finite-size scaling hypothesis for Sn, i-e. assum- 
ing that SniCi N) takes a homogeneous scale-invariant 
form with respect to its arguments in the vicinity of the 
transition point, ^c, when both N and A are sufficiently 
large. In other words, for any positive number A > 0, 

5„(AC, A-'-A, A-^iV) = X'^SniC, A, N), (88) 

where a, b and c are a new set of exponents associated 
with A and N respectively. The above relation implies 
that when the reduced Manning parameter, C, is rescaled 
with a factor A, the size parameters, N and A, can be 
rescaled such that the order parameter remains invariant 
up to a scaling prefactor. Finite-size scale-invariance is 
a common feature in critical phase transitions |49l It^ 
and provides an accurate tool to estimate critical expo- 
nents in numerical simulations |75l . The exponents 
in Eq. H88|) can be calculated directly from MC simu- 
lations. These exponents are in fact related to and give 
the values of the desired critical exponents 7, v and X; 
as will be shown below. Note that the exponents may 
in general depend on n (the index of Sn), the coupling 
parameter, S, or the space dimensionality, which are not 
explicitly incorporated in the proposed scaling hypothe- 
sis, but their influence will be examined later. 

Given Eq. (|88|) , the following relations are obtained by 
suitably choosing A. For A — N^^^, one flnds 

5„(C, A, N) = N-^'^CniCN^'", AiV-''/"), (89) 

where C„(m, v) is the scaling function corresponding to a 
system with both finite N and A. The above expression 
is useful for a system with finite N in the limit A — > cxj. 
Thus assuming that Cn{u,v) exists for v = AN^''^'^ — > 
00, the relation reduces to 



^„(C,oo,iV) =iv-''/W„(Civi/=), 



(90) 
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where the scaling function JVn{u) = C„(u, oo). The crit- 
ical exponent v follows by considering this relation right 
at the critical point, ( — 0, i.e. 



^„(0,(X3,iV)=AA„(0)7V-^ 
where v is obtained as 



a 

u — —. 

c 



(91) 



(92) 



On the other hand, we assume that in the vicinity of 
(and above) the critical point (i.e. for small but finite 
C), Sn{<^,oo,N) is only a finite function of the reduced 
Manning parameter C, when the limit — > oo is taken. 
Hence the scaling function Mn{u) is required to behave 
as A/'„(u) ^ It" for m — > oo, which yields 

5„(C,oo,(^)^C^ (93) 

where the critical exponent associated with C reads 

X = a. (94) 

To determine the critical exponent associated with A 
in terms of the exponents {a, 6, c}, we need to consider 
Eq. (jHHl) for A = Ai/''. We thus have 



5„(C, A,7V) = A-'^/''C'„(CAl/^iVA-=/^), 



(95) 



where C'n{u,v) is a new scaling fimction. This relation 
is useful for a system with finite A in the limit N oo, 
where assuming that C'n{u,v) exists, we obtain 



5„(C,A,oo) = A-'^/^I?„(CAi/'') 



(96) 



with a new scaling function 'Dn{u) = C„(u,cx)). The 
critical exponent 7 follows by considering this relation 
right at the critical point, C, —0^ that yields 



5„(0,A,c5o) =P„(0)A- 



where 7 reads 



(97) 



(98) 



Therefore, we have a complete set of relations (|92ll . (|94f) 
and H98() from which the critical exponents 7, v and x 
may be obtained using the exponents a, h and c. 

Equation (|96|) compares directly with the mean-field 
result, Eq. where we showed that 7pb = 2 and 

Xpb ~ 2. Note also that the exponent v is not defined 
within the mean-field theory. 



D. Critical exponents: the CCT universality class 

1. The exponents x o,f^d u 

In order to verify the generalized finite-size scaling hy- 
pothesis (|88|l and estimate the critical exponents numer- 
ically, we shall adopt the standard data-collapse scheme 
used widely in literature |77| . 
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FIG. 11: Rescaled order parameter, N'^^'^Si, as a function 
of the rescaled reduced Manning parameter, (^N^^'^, in the 
vicinity of the critical point, = 1.0, and for small and large 
coupling parameters H = 0.1 (main set) and S = 10'^ (inset). 
Symbols show data for various number of particles A'^ = 50 
(triangle-downs), 70 (circles), 75 (squares), 100 (diamonds), 
200 (cross symbols), 300 (triangle-ups) , and fixed A — 300. 
In these plots, the exponents are chosen here as a/c — 2/3 
and 1/c = 1/3. Error-bars are smaller than the symbol size. 



We begin with the exponents x ^-nd v that can be cal- 
culated using Eq. (|89|l . which involves a scaling function, 
Cn{u,v), of two arguments u — QN^^'^ and v — AN~''^'^. 
In our simulations, A'' ranges from 25 up to 300 and A 
ranges from 50 up to 300; thus assuming that the expo- 
nent b/c is small, which will be verified later on, we deal 
with a typically large value for v ^ 10 ~ 10^. There- 
fore the limiting relation H90|l is approximately valid and 
yields 



(99) 



Now if the data for Sn are plotted as function of ^ = 
1 — ^c/C for various N (but at fixed sufficiently large 
A), equation H99|) predicts that by rescaling the reduced 
Manning parameter C by the factor of N^^'^ and the order 
parameter by the factor 7V°/'^, all data should collapse 
onto a single curve. Numerically, this procedure allows to 
determine the exponents a/c and 1/c in such a way that 
the best data collapse is achieved. We show the results in 
Figure 1111 for various N (symbols) and for the coupling 
parameter S = 0.1 (main set). The collapse of the data 
onto each other is indeed achieved within the numerical 
error-bars for the exponents in the range 1/c ~ l/3±0.05 
and a/c ~ 2/3 ± 0.1. This yields the critical exponents 
and X from Eqs. (|^ and as 



i^~2/3±0.1, 
2.0 ±0.4, 



(100) 
(101) 



where the errors are estimated using the standard error 
propagation methods. The value obtained for x agrees 
with the mean-field result, Eq. (|44|l . 
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FIG. 12: Main set: Order parameter Si as a function of the 
lateral extension parameter, A, for increasing number of coun- 
terions from A*' = 70 up to 300 as indicated on the graph 
(Manning parameter here is ^ = 1.0 and the coupling param- 
eter S = 0.1). Dashed line shows the PB power-law, Eq. (1451 . 
Inset: Rescaled order parameters Sn = (j^ ") as a function 
of n for Manning parameters close to the critical point (from 
top): ^ = 1.03 (filled circles) and ^ = 1.01 (filled diamond). 
Here the coupling parameter is H = 0.1, the number of coun- 
terions N = 100, and A = 300. 
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FIG. 13: Rescaled order parameter, A°'^''Si, as a function 
of the rescaled reduced Manning parameter, ^A^'"', in the 
vicinity of the critical point, ~ 1-0, and for small and large 
coupling parameters H = 0.1 (main set) and H = 10^ (inset). 
Symbols show data for various lateral extension parameters 
A = 100 (circles), 150 (cross symbols), 200 (diamonds), 250 
(plus symbols) and 300 (triangle-ups) . Number of counterions 
is fixed {N = 200 for the coupling parameter H = 0.1 and 
TV = 100 for H = 10^), and the exponents are here chosen as 
a/6 = 2.0 and 1/6 = 1.0. 



In order to check whether the exponents vary with the 
electrostatic coupling parameter, S, we repeat this pro- 
cedure for a wide of range of values for S. We find the 
same values for the exponents for coupling parameters up 
to S = 10^. For comparison, the results for S = 10'^ are 
shown in the inset of Figure ^] where the data collapse 
is demonstrated for l/c — 1/3 and a/c= 2/3. 



2. The exponent 7 

Given the exponents a and c calculated above and mak- 
ing use of the finite-size scaling relation , we can esti- 
mate the exponent 6, and thereby the scaling exponent 7, 
associated with the lateral extension parameter. In this 
case, however, the second argument v — N/Sr'^/^ in the 
scaling function C'n{u,v) defined in Eq. may not 

be considered as large within our simulations (since as 
shown below the ratio c/b is large). But it turns out that 
the dependence of C'niu,v) on v is quite weak such that 
the finite-size scaling relation H96|) is approximately valid 
and can thus give the desired exponent. To examine this 
latter property of C'„(u, v), we consider the relation H95() 
right at the critical Manning parameter ((^ = 0), i.e. 

S'„(0, A, N) = A-^/^CniO, NA-"/''). (102) 

In Figure [T^ <S'n(0, A, N) is plotted as a function of A in 
a log-log plot for increasing N from 70 up to 300 and for 
S = 0.1. As clearly seen, the order parameter varies quite 
weakly with the number of particles, and the variations 



are already within the error-bars (equal to symbol size) 
for N > 100. 

Thus multiplying both sides of Eq. H9t)|) with A"/'', we 
have 

A^/'Sn ^ 2?„(CAi/'), (103) 

in which the exponent a previously determined as a = 
X = 2.0 ± 0.4. We thus plot the order parameter Sn as 
a function of ^ for various A (but at fixed sufficiently 
large N) , and rescale both Sn and C values with the scal- 
ing factors A"/'' and A^/'' respectively; the exponent b 
is chosen in such a way that the best data collapse is 
obtained within the error-bars. The result is shown in 
Figure 1121 for A'^/^'S'i as a function of CA^/*", where the 
coupling parameter is chosen as S = 0.1. The collapse of 
the data onto each other is obtained only for the expo- 
nent 1/fe in the range l/&~ 1.0±0.2, yielding the critical 
exponent 7 from Eq. H98|) as 

7 -2.0 ±0.6, (104) 

which agrees with the mean-field exponent, Eq. H47|l . 
We find the same value for 7 by repeating the above 
procedure for larger coupling parameters. For instance, 
the results for S = 10^ are shown in the inset of Figure 
1131 where we have chosen l/b — 1.0. 

Note that the estimated values of b and c show that 
the ratio b/c is as small as 0.3, which is consistent with 
the assumption made in using the asymptotic forms H90() 
and H9t)|l in the foregoing data-collapse procedure. 

As a main result, our numerical data confirm the ex- 
istence of characteristic scaling relations associated with 
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the counterion-condensation transition in 3D and show 
that the values of the critical exponents are universal, 
i.e. independent of the coupling parameter, S, and agree 
with the mean-field universality class. 

Also, in agreement with mean-field results, the expo- 
nents are found to be independent of n, the index of the 
order parameters Sn = (1/f"). In fact, we find that the 
higher-order moments are related to the first-order mo- 
ment, 5*1, via 



n 



(105) 



in the vicinity of the critical point, which indicates that 
nSn is independent of n, as demonstrated in the inset of 
Figure ^1 (compare with the mean-field relation H41|l ). 



VII. COUNTERION-CONDENSATION 
TRANSITION (CCT) IN TWO DIMENSIONS 

In this section, we shall investigate the role of space 
dimensionality in the asymptotic behavior of counteri- 
ons at a charged cylindrical boundary, by considering 
a 2D counterion-cylinder model. As a typical trend in 
bulk critical phenomena, the effects of fluctuations near 
the critical point are known to grow with diminishing di- 
mension 49], and cause large deviations from mean-field 
theory. It is therefore interesting to study the CCT in a 
lower spatial dimension. 



A. The two-dimensional model 

In 2D, we use a primitive cell model similar to the 
3D model described in Section III Al It consists of a 2D 
central charged cylinder (central "disk" ) of radius R con- 
fined co-axially and together with its neutralizing point- 
like counterions in an outer cylinder (outer "ring") of 
radius D. In order to construct the two-dimensional in- 
teraction Hamiltonian, we use the fact that the Coulomb 
interaction between two elementary charges in 2D (the 
2D Green's function) is of the form 



1'2d(x) 



Inlxl. 



(106) 



This follows directly from the solution of the 2D Poisson 
equation for a point charge, that is 



(107) 



The configurational Hamiltonian of the 2D system may 
thus be written as 



n 



N 



N . s 



R 



(108) 



with Xi = {ri,4>i) being the position vector of the i-th 
counterion (in polar coordinates), and Ac and being 



dimensionless charges of the counterions and the cylinder 
respectively. The first term gives the counterion-cylinder 
attraction and the second term gives mutual repulsions 
between counterions. Clearly, the present 2D model is 
equivalent to a 3D system comprising an infinitely-long 
central cylinder (of radius R) in the presence of mobile 
parallel lines of opposite charge as "counterions" , which 
may be applicable to a system of oriented cationic and 
anionic polymers [s^ ■ Using this 3D analogy, the prefac- 
tors Xr and Ac may be related to the linear charge density 
of the cylinder and counterion lines respectively. 

Taking the logarithmic interaction H106|) will also en- 
sure that the general form of the field-theoretic represen- 
tation for the system remains the same as in the 3D case 
|78j| . and in particular, the mean-field Poisson-Boltzmann 
theory, which follows from a saddle-point analysis, is rep- 
resented exactly by the same equations and results as 
discussed in Section Hill 



B. Rescaled representation 

In analogy with the 3D system, we shall refer to the 
dimensionless prefactor of the counterion-cylinder inter- 
action in Eq. (|1U8|I as the Manning parameter, that is 



C — AcAj./2. 



(109) 



Also the prefactor of the countcrion-counterion interac- 
tion is defined as the coupling parameter 



S = A^/2. 



(110) 



These definitions can be justified systematically when the 
Hamiltonian of the system is mapped to an effective field 
theory, where S and ^ formally appear in the same role 
as in 3D [t^ ■ We shall conventionally rescale the spatial 
coordinates as x = x//i2D using the length scale = 
R/^, which is the 2D analogue of Eq. 10. 
The Hamiltonian in rescaled units reads 



n 



N 



kBT 



2e 



N 



R 



2-5: In 



R 



(111) 



The electroneutrality condition implies A^ = NX^, where 
N is the number of counterions in the system. This re- 
lation may also be written as 



(112) 



Thus an important consequence of electroneutrality in 
2D is that the coupling parameter and the Manning pa- 
rameter are related only via the number of counterions. 
In particular, in the thermodynamic limit N — )■ oo, the 
coupling parameter tends to zero, S ^ 0, suggesting that 
the mean-field prediction should become exact! 

We use a similar simulation method as devised for the 
3D system using the transformed coordinates (y, 0) with 
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y = ln{f/R) being the logarithmic radial distance of par- 
ticles from the central cylinder. As explained in Section 
Ivl this transformation leads to the centrifugal sampling 
method appropriate for equilibration of systems with 
large lateral extension parameter A = \n{D/R) ^ 1, 
where the critical behavior associated with the CCT 
emerges. The 2D partition function thus reads 



R 



2N 



-■N 



m 



V 



N 



d0i dyi 



exp 



n 



N 



(113) 



where < ?/ < A and the transformed Hamiltonian, 

N 



H 



N 



(2C-2)^y.-2S^ln 

i=l (ij) 



R 



(114) 



The minimal set of dimensionless parameters in 2D is 
given by the Manning parameter, ^, total number of 
counterions, N , and the lateral extension parameter A. 
The range of simulation parameters and other details are 
consistent with those given in Section IV Bl 
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FIG. 14: Order parameter Si — (1/r) as a function of Man- 
ning parameter, ^, for the 2D counterion-cylinder system. 
Symbols show simulation data for different number of par- 
ticles as indicated on the graph. The mean-field PB (solid 
curve) and the strong-coupling (dashed curve) predictions are 
obtained from Eqs. 1401 and 1671 respectively, which are 
valid in 2D as well. The lateral extension parameter here 
is A = 300. Thin dashed curves are guides to the eyes. 



VIII. SIMULATION RESULTS IN 2D 
A. The order parameters 

We consider the same set of order parameters Sn = 
(1/f") as defined in Eq. H38|) to characterize the CCT in 
2D. They can be measured in the simulations as 

i=l 

for n > 0, where the bar sign denotes the MC time aver- 
age after proper equilibration of the system. Of partic- 
ular interest is the behavior of Sn as a function of Man- 
ning parameter, ^, which identifies the two regimes of 
complete de-condensation (with vanishing Sn) and par- 
tial condensation (with > 0) as A ^ oo. Unlike in 
3D, where S can be varied as an independent parameter, 
various coupling regimes in the 2D system are spanned 
by changing the number of particles, N, for a given ^ (see 
Eq. mUl ). 

The 2D simulation results for the order parameter Si 
are shown in Figure [HI for increasing number of particles 

= 1, 2, 3, 5, 10 and 100 (symbols) and for a large lateral 
extension parameter A = 300. As seen for the smallest 
number of counterions, A^ 1, the data trivially fol- 
low the strong-coupling prediction, Eq. H67() . shown by 
the dashed curve (Section II VI) . For increasing A^, Si 
decreases and for sufficiently large values, the data con- 
verge to the mean- field PB prediction, Eq. H40|l . shown 
by the solid curve. This in fact occurs for the whole 
range of Manning parameters and thus confirms the trend 
predicted from the 2D electroneutrality condition (|112|l . 
Accordingly, scaling analysis of the order parameters for 



large N gives identical results for the critical exponents 
as in 3D fSections lVI Cl and lVI D|l and thus in agreement 
with the mean-field theory, which we shall not discuss 
here any further. The result that the mean-field the- 
ory for the counterion-cylinder system is exact in 2D for 
A^ — > cxD is in striking contrast with the typical trend in 
bulk phase transitions "4^ , and also with the situation in 
3D, where the strong-coupling effects become important 
in the condensation phase (^ > 1) for growing S (Section 
IVIA4^ . 

The order-parameter data in Figure 1141 on the other 
hand, reveal a peculiar set of cusp-like singularities, 
that are quite pronounced for small number of particles. 
These points become strictly singular in the limit A cx3 
and represent the Manning parameters at which individ- 
ual counterions successively condense (or de-condense). 
We will demonstrate this point using an analytical ap- 
proach in Section IVIII CI (A similar singular behavior 
is also found in 3D for small A^, but the behavior in 3D 
appears to be more complex and will not be considered 
in this paper). 



B. Energy and heat capacity 

The singularities at small particle number, A^, appear 
also in the internal energy and the heat capacity. In 
Figures El and ^1 we plot the rescaled energy, E — 
En / [Nk^T), and excess heat capacity, C = Cm /{Nk^), 
obtained from the simulations using Eqs. (|78|l and H79() 
and the 2D Hamiltonian (|111|) . as a function of ^ and 
for A^ = 1,2,3,4 and 5. As seen, the energy shows a 
sawtooth-like structure for increasing ^ consisting of wide 
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FIG. 15: The rescaled (internal) energy of the 2D counterion-cyhnder system, E — En /{NksT), as a function of Manning 
parameter, ^, for different number of particles as indicated on the graph. Singular regions represent successive condensation 
(de-condensation) of counterions-see Section [Vlll C II The lateral extension parameter for these data is A = 300. The dashed 
curves are the analytical results given by Eq. 112411 . 
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FIG. 16: The rescaled excess heat capacity of the 2D system, C = C n / {NkB), as a function of Manning parameter, ^, for 
different number of particles, A'^, as indicated on the graph (for clarity, the data are here multiplied by A''). The peaks represent 
successive condensation (de-condensation) of counterions-see Section rVIII C II For these data A = 300. The dashed curves are 
the analytical results given by Eq. 11125^ . 



regular regions, in which the energy almost linearly in- 
creases, and narrow singular regions, where the energy 
rapidly drops. Recalling the thermodynamic relation, 




it follows that the excess heat capacity vanishes in the 
regular regions, but develops highly localized peaks in the 
singular regions, as also seen from the simulation data in 
Figure El 



C. Condensation singularities in 2D: an analytical 
approach 



In what follows, we present an approximate (asymp- 
totic) analysis of the 2D partition, which elucidates the 
physical mechanism behind the singular behavior in 2D. 
The rigorous analysis of the 2D problem is still missing 
and more systematic approximations have been devel- 
oped recently [t^- 
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1. The partition function 

Suppose that the Manning parameter is such that N — 
m counterions are firmly bound to the central cylinder 
(disk) , while m counterions have de-condensed to infinity, 
where m — 1, . . . ,N. Using the 2D Hamiltonian 
the partition function can exactly be written as 



Z 



N 



N 



n 



exp 



{-TS^}'<n4?(n7) 



in actual units, where T-LN-m represents interactions 
among condensed counterions (labeled by i — m + 
1, ... ,7V), and 



r(0 

-■N 



dx; exp < — 2^ In 



2C 

N 



N 



R 



(118) 

is the contribution from individual de-condensed counte- 
rions (labeled by / = 1, . . . , m). Assuming that the de- 
condensed counterions are de-correlated form each other 
and also from the condensed counterions as they diffuse 
to infinity for A — > oo (i.e. using jx; — x;| ~ r;), Z 
approximately factorizes as 



(0 

N 



^ ^27r/ndr,exp^-2ein(^')-f| E In ^ ""^ 



ril) 



2t:R- 



exp [2(1 -^//iV) A] - 1 
2(1 -eV^) ■ 



R 
(119) 



In the limit A — s- cx), Z]^ diverges for Manning parame- 
ters 

N 

T 



(120) 



which indicates de-condensation of the l-th counterion 
from the charged cylinder (see Section FlI D|) . Repeating 
the above argument for various number of de-condensed 
counterions, one finds a set of singular Manning param- 
eters. 



N 

T 



for I 



(121) 



at which individual counterions de-condense from the 
charged cylinder. These singular points coincide with 
the values obtained from our simulations based on the 
full partition function H113|) for very large A (see Figures 
I14ll6l and Table for the numerical values). 



2. Energy and heat capacity 

In general the partition function H117|l can exactly be 
written as 



JV 



(122) 



1=1 



where Z^ is defined in Eq. H118|l . For A » 1, the dom- 
inant contribution to the internal energy and the heat 
capacity comes from de-condensed counterions. Thus, in 
order to derive analytical expressions for energy and heat 
capacity on the leading-order for large A, we shall use Eq. 
(|122|l together with the approximate expression (|119|l for 

Z^\ Hence, we obtain the leading-order contribution to 
the free energy J^N/{kBT) = — IuZat as 



1=1 



" g2(i-e/^f)A _ ^ 



NkBT 



(l-^/Cf) 



(123) 



for A ^ 1, and thereby the (rescaled) internal en- 
ergy E — £_dJ^/d£, and the rescaled heat capacity C = 
—(^d'^T/d^'^ are obtained as 
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(124) 



C: 
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sinh^[2(l-e/ef)A] 



(125) 

The above expressions are shown in Figures 1151 and 1161 
for A — 300 and for various number of particles (dashed 
curves), which as seen closely reproduce the behavior ob- 
tained in the simulations (symbols). 

Note that as an individual counterion de-condenses at 
^f, the internal energy suddenly jumps, since the de- 
condensing counterion gains a large amount of energy 
due to its logarithmic interaction with the central cylin- 
der. The regular regions (between two successive jumps) 
in the energy curve are dominated by de-condensed coun- 
terions and thus exhibit linear scaling with A = \ii(D/R). 
The asymptotic form of the energy in these regions for 
A ^ oo is obtained from Eq. (|124|l as 



E + 
lim — — TT 

A->oo A 



1) 



C for ef+i < C < 



(126) 



The singular part of the energy corresponds to a nar- 
row region around each : which (except for the upper- 
most singularity) is bounded between a local minimum 
(slightly above ) and a local maximum (slightly below 
it)- The approximate locations of these extrema are ob- 
tained as 



1 + 



1 



and 
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1 



(127) 

using Eq. (|124|l . and for large A. The energy jump, 6Ei, 
upon de-condensation of a counterions at ^ = is then 
given by 



dEi = Eurn 



2A 



(128) 
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Note that this value can also be obtained directly from 
Eq. H126|l . For A ^ oo but at fixed and finite N, the 
energy curve tends to a sharp sawtooth-like form as both 
f""" and ^["^'^ tend to producing TV strictly step-like 
singular points, at which the limiting energy jump is 
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(129) 



The heat capacity expression (|125|l . on the other hand, 
exhibits N isolated peaks for A ^ 1. The heat capacity 
at diverges as C{£, = ^f) ~ A^/3 with increasing A, 
giving rise to N limiting algebraic divergencies as 



N 

lim C^V 

1 = 1 



(130) 



D. Critical point and the continuum limit 

The lower-most singularity located at ^ = is asso- 
ciated with the de-condensation of the "last" counterion 
from the charged cylinder. As shown above, this singu- 
larity occurs at unity (^^ = 1) when A — s- oo and is thus 
independent of the number of counterions. It therefore 
gives the exact location of the 2D critical point as = 1 
when the continuum (thermodynamic) limit iV — > cxd is 
approached, which coincides with the mean-field predic- 
tion. Note that in analogy with the method used in Sec- 
tion IVI Bl can also be derived from the asymptotic 
value of the energy maximum location, Eq. H127|) . for 
I = N, when A and N both tend to infinity. 

Equations p28|l - H130|l represent the asymptotic results 
when the system size increases to infinity but the num- 
ber of particles, A'^, is finite. In the converse limiting 
case, i.e. when A is large and fixed but N increases to 
infinity (continuum limit), all singularities smoothen ex- 
cept for the one, which represents the critical point. The 
limiting energy curve for N oo may be obtained as 
follows. First note that the width of the energy jump 
around each singularity tends to zero as indicated by Eq. 
(|128|l . Secondly, the spacing between singular points 
(and thus the width of regular regions for C > 1) tends to 
zero (as ~ ^/N) with increasing N. Therefore, the energy 
at a given Manning parameter ^ between two successive 
singularities, ^^^-^ < ^ < ^j", is approximately given by 
E ~ E{(^ = £^i), where the right hand side is obtained 
from Eq. ifT^ a.s E{£, ^ ) = A/^f . This imphes that 
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for ^ > 1 and sufficiently large A. For small Manning 
parameter ^ < 1, there are no singularities and from Eq. 
(|124|l . we obtain 
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TABLE I: Numerical values of the location of the singular- 
ities, Eq. 112H . for the 2D counterion-cylinder system for 
several number of particles (compare Figures Il4ll6ll . 



for large A. These limiting results can also be obtained 
using Eq. ifT^ . 

The energy curve in the continuum limit therefore coin- 
cides with the universal form obtained within the mean- 
field theory in Section UTIEll (see Eq. (EU). The heat 
capacity in this limit follows from Eq. (|81|l . and exhibits 
a universal jump at = 1 in agreement with Eq. (|62|l . 



E. The condensed fraction 

The preceding results enable us to calculate the limit- 
ing condensed fraction of counterions, a(^), as well when 
A — > oo and N —^ oo. For a given Manning parameter, ^, 
and number of particles, N, the condensed fraction Q;Ar(^) 
is given by the number of singularities located below ^, 
i.e. 



ajv(0 



I 



for er+i < e < cf 



(133) 



This fraction is trivially zero for < = 1 as A oo. 
Using Eqs. (|121|l and (|133|l . we obtain the condition 

- ^ < 1 - ^ < "^(0, (134) 
which in the limit of infinite number of counterions yields 



«(e)= lim ajv(e) = l-7- (135) 

This is nothing but the mean-field or Manning condensed 
fraction, au = 1 — 1/^ (Section IVI A2|) . The finitc-ion- 
number correction to this limiting value follows from Eq. 
(Pl| as 



(136) 



IX. CONCLUSION AND DISCUSSION 

In this paper, we present an extensive numerical anal- 
ysis of the critical behavior of counterions at a charged 
cylinder (counterion-condensation transition) in both 
two and three spatial dimensions. Analytical results for 
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the critical behavior are also derived using the asymp- 
totic theories of mean field (Poisson-Boltzmann) and 
strong coupling. The counterion-condensation transition 
(CCT) is regulated by the dimensionless Manning pa- 
rameter (rescaled inverse temperature), ^ = q^BT, and 
occurs at a critical threshold ^c, below which counterions 
completely unbind (de-condense) to infinity, but above 
^c, a finite fraction of counterions binds (or condenses) 
in the vicinity of the charged cylinder. Since the CCT 
criticality emerges asymptotically in the limit of infinite 
system size and infinitely many particles, we have em- 
ployed Monte-Carlo (MC) simulation of a periodic cylin- 
drical cell model in logarithmic radial coordinates, which 
gives rise to a powerful (centrifugal) sampling method 
for extremely large lateral system sizes within reason- 
able equilibration times. This constitutes the key part 
of the present numerical investigation, since the critical 
and universal aspects of the CCT within the cell model 
can only be captured for large logarithmic system size 
A — \n{D/R) ^ 1 (with D and R being the outer bound- 
ary and the charged cylinder radii respectively). 

As main results, we have determined the precise loca- 
tion of the critical Manning parameter, ^c, the scaling 
universality class of the CCT and the singular behav- 
ior of energy and heat capacity on a systematic level 
and without suppressing inter-particle correlations. As 
shown both the mean internal energy and the excess heat 
capacity become singular at the critical point. The ex- 
cess heat capacity, which vanishes in the de-condensation 
phase, shows a universal discontinuity (jump) at the crit- 
ical point indicating that the CCT is a second-order tran- 
sition, as also suggested in a recent mean- field study 's^ . 
In a finite system, these singularities appear in the form 
of a pronounced peak, the asymptotic behavior of which 
is used to determine the critical Manning parameter, ^c, 
in the simulations. On the other hand, the critical ex- 
ponents associated with the CCT are obtained using a 
combined finite-ion-number, and finite-size. A, anal- 
ysis of the order parameters 5„ = {r^") (with f = r/fi 
being the radial distance from the cylinder axis in units of 
the Gouy-Chapman length, /x). These order parameters 
represent the mean inverse localization length of counte- 
rions. For ^ < and in an infinitely large system, Sn 
vanishes, but takes a finite value above ^c, which exhibits 
the scaling relation Sn ^ C^, where C = 1 — ^c/C is the 
reduced Manning parameter (the reduced temperature) 
and the exponent x is determined as x = 2.0 ±0.4 within 
our MC simulations. In finite systems, Sn does not van- 
ish at and displays a power-law decay with increasing 
size parameters, A and N, as S'n(C = Cc) = A~'^ (when 
number of particles, iV, is fixed) and S'„(^ = ^c) = 
(when lateral extension parameter. A, is fixed), where 
the critical exponents are determined as 7 = 2.0 ± 0.6 
and v = 2/3 ±0.1. 

The critical exponents are demonstrated to be univer- 
sal, i.e. independent of the coupling strength, S (varied 
over several decades 0.1 < S < 10^), and agree with 
the values we obtained from the mean-field PB theory 



as xpb = 2.0 and 7pb = 2.0 (note that the exponent v 
is not defined within mean-field theory). Interestingly, 
the critical Manning parameter both in 3D and 2D is 
also found to be universal and given by the mean-field 
value = 1- Therefore, in striking contrast with the 
typical situation in bulk critical phenomena, where de- 
viations from mean-field theory due to fluctuations grow 
with decreasing space dimensionality, the CCT critical- 
ity is found to be described by the mean-field universal- 
ity class both in 3D and 2D. Surprisingly, the mean-field 
theory is found to be exact in 2D in the whole range of 
Manning parameters (or temperatures) when N ^ 00. 
In 3D, however, correlation effects in the condensation 
phase {S, > S,c) lead to strong deviations from mean- 
field theory and support the strong-coupling predictions 
characterized by an excessive accumulation of counteri- 
ons near the cylinder surface in agreement with previous 
numerical studies 0, l40l l6^ 1651 l66l | and experiments 
[73] . An important result is that the large-distance form 
of the density profile remains unaffected by these cor- 
relations and thereby a universal condensed fraction is 
obtained when the inflection-point criterion is applied. 

To our knowledge, rigorous analytical derivation of the 
critical Manning parameter or the scaling exponents of 
the CCT in 3D is not yet available and our simulations 
provide the first numerical results for the universal and 
critical features of this transition in the large-system-size 
limit In 2D, we have shown that the simulation 

results can be understood using an approximate analyt- 
ical theory. The present predictions for order parame- 
ters and thermodynamic quantities can be examined in 
experiments. In particular, the order parameters may 
be obtained from the distribution of counterions around 
charged polymers, which has been directly measured us- 
ing anomalous scattering techniques |80l |. 

In this study, we have made use of an idealized (primi- 
tive) cell model [sollsills^ in order to bring out the main 
universal aspects of the CCT. It is interesting to examine 
possible effects due to more specific factors that exist in 
realistic situations, namely, the discrete charge pattern 
of polymers H IM IH IsaJSJel^^ flexibility 
and finite contour lengthlunillllE IH IH IH Hi 
I41I Isil Is^ l as well as the infiuence of non-uniform dielec- 
tric boundaries |4ll | on the critical behavior. However, 
the present results already indicate that short-range ef- 
fects such as electrostatic correlations do not affect the 
properties of the system near the critical region (^ ^c), 
since most of counterions are de-condensed and the crit- 
ical behavior is predominately determined by long-range 
features. Future studies will be useful to address the 
possible influence of short-range specific effects on the 
critical behavior [83j |. 

In this work, we have not investigated the role of addi- 
tional salt and co-ions, which lead to screening of electro- 
static interactions (see, e.g., Refs. 0,111; EM E3)- It is 
known that the Debye screening length, r^, plays the role 
of an upper bound cut-off (similar to the outer boundary 
D in the cell model): the CCT occurs for vanishing salt 
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concentration, i.e. when ln(rs/i?) ^ oo 0, 0, 0, l43l |. 

Thus we can expect similar asymptotic behaviors arise 
near = 1 and within a similar model as used here 
when the vanishing-salt limit is approached [s^ . 
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APPENDIX A: SINGULARITIES ASSOCIATED 
WITH THE ONSAGER INSTABILITY 

The rescaled energy, E — EN/{NkBT), and excess 
heat capacity, C — Cn / (Nk-B), associated with the On- 
sager instability can be calculated from Eq. (|17|l . The 
results coincide with those obtained in Section IVlII GI bv 
choosing = 1 (see also Figures ITHland lT^ . In brief, the 
approximate location of the energy peak is obtained 
as 



1 



1 



(Al) 



for finite large A, which shows a different asymptotic 
convergency (from below) to the critical value = 1 as 
compared with the mean- field result, Eq. (jH2J- The heat 
capacity peak is located above the critical point at 



1 



5 

A2- 



(A2) 



In the limit A ^ oo, the heat capacity diverges at the 
critical point from above and below as 



(A3) 



where ^ = 1 - 1/^. The left tail of energy (for ^ < ) 
goes to infinity linearly with A as E ~ 2^A-|-^/(^ — 1), 
but its right tail shows a power-law divergency as 



E^C 



(A4) 



Note that these behavior are distinctly different from 
those obtained in the simulations with many particles 
fSection lVI B|) and within the mean-field theory (Section 
UTTTT^ . 



APPENDIX B: RESCALED PB EQUATION 

The Poisson-Boltzmann equation for mean electro- 
static potential ^cicc in actual units reads 



V^V'elec 



(T(x)e qep{'x) 



(Bl) 



where the density profile of counterions is given by 

p(x) = po f^(x) e-«^'3'^-"^ (B2) 

with po being a normalization prefactor and ri(x) spec- 
ifying the volume accessible to counterions (Eq. 120() '). 
The rescaled PB equation l|18|) is recovered using above 
equations and the dimensionless parameters as x = x//i, 
ct(x) = P(t(x)/cts, /5(x) = p(x)/(27r£BO's )i and also 
i) = qeV'cioc/(fcBT), where p = R/i = l/(27r(j^Ba's) is 
the Gouy-Chapman length (Eq. 10). The prefactor 
po is related to k in Eq. H18(l through k — Kfi where 
= 47rq^£BPo- One can show that the density profile in 
Eq. ljB2|l is also mapped to the rescaled density H23() in 
the text. 

Normalization of the density (to the total number of 
counterions) in actual units reads / dx p(yi) = N, and 
in rescaled units, we have J dxp(x) = 27rSA^, which is 
equivalent to Eq. (|21|l in the text, when Eqs. Hll|) and 
(1^ are used. 



APPENDIX C: ASYMPTOTIC RESULTS 
WITHIN THE PB THEORY 

1. Limiting behavior of (3 for large A 

A general discussion has been given by Fuoss et al. [HJ 
for the overall behavior of /3 using Eq. (j26(l . Here, we 
first review some of their results as quoted in the text, 
and then derive the finite-size scaling relations for f3 near 
the PB critical point (^^^ = 1) as used in Section UlI Gl 



a. Small Manning parameter ^ < Aaf 

The integration constant /3 vanishes at ^ = Aaf and 
tends to unity, /3 ^ 1~, for small ^ 0"*", as it can 
be checked easily from Eq. H26|l (we arbitrarily choose 
/3 > 0) ^El\. Further inspection shows that in this regime, 
/3 ^ (^ — 1)^ when A oo Hence for A ^ 1, one 
can propose the following form 



f3^ ^i^-lfil-x), 



(Gl) 



where x is a small function of ^ and A. To determine x, 
we may rearrange the first equation in (|26|l as 



/^A^ilni^-iln(i-l4±^. 

2 l + (3 2 (e-l)-/3 

and use this together with Eq. (jGip to obtain 



4C 



=2(?-l)A 



2-C 

which reproduces Eq. (|31|l in the text. 



(G2) 



(G3) 



27 



b. Large Manning parameter ^ > Aap .• 

For large Manning parameter ^ > Aaf, Eq. ((JSJ, 
tends to a finite upper bound Poo = t^/^, when ^ — > oo 
[Slf . Thus for A ^ oo , /3 vanishes for the whole range of 
Manning parameters > Aaf — 1 as used frequently in 
the text (see e.g. in Eq. lH^ I. 



c. Finite-size scaling for (5 close to 

Of particular importance in our analysis is the behavior 
of P close to = 1. (Since always Aaf < 1, we restrict 
the discussion only to the regime ^ > Aaf)- Analysis 
of Eq. I|26(l shows that for sufficiently large A, we have 
/3 ~ 7r/(2A) right at the critical point = 1. We may 
then perform a Taylor expansion around to obtain 
the approximate form of (3 for small C = 1 — ^J'^/^. We 
find 



/3(C,A) = - + 



2 , 8A 

-c- 

TT TT 



3 + o{e), (C4) 



which remains valid for ^A < 7r^/4. This relation clearly 
indicates a scale- invariant form for /3 when A ^ 1. Com- 
paring this with Eq. (|50|1 . we find the approximate form 
of the scaling function B{u) as 



(C5) 



where u = CA < 7r^/4. In particular, we have B{u) 
tt/2 as u 0. 

The asymptotic behavior of B(u) for u "> cxd (or equiv- 
alently A — > oo for finite C) can be obtained using a 
different series expansion, since in this limit /? becomes 
singular at — 1 and the above expansion breaks 
down. This is because P is always singular (with an in- 
finite first derivative) at Aaf which tends to the critical 
Manning parameter [5l| . We thus perform an expansion 
around ^ = Aaf, which yields P ~ a/S^/A. This gives 
the asymptotic form of the scaling function as 



Biu) 



3u 



oo. 



(C6) 



2. The PB cumulative density profile 



The PB cumulative density of counterions, rtpB(y), 
is defined through Eq. (|73|l . It can be easily shown 
that npB{y) is a monotonically increasing function of 
y = ln(r/i?), i.e. dnpe/dy > |23|. It is there- 
fore bounded by its boundary values npB(O) = and 
npB(A) = N (see Figure EJ. 

For ^ > Aaf, npsijl) has an inflection point the lo- 
cation of which, — ln(r,/i?), follows from equation 
d^Tipe/dy^ = as = tan~^[(^ - l )//3]/ /3. It is easy to 
check (using the results in Appendix IC l|l that ~ A/2 
for large ^ and that y» for ^ — > 1 (y* becomes nega- 
tive for smaller values). For ^ < Aaf, on the other hand. 



the cumulative density, ?t.pb(j/), vanishes for y < A as 
A — > cxD, as can be checked by inserting the approximate 
expression (ICTl for P into Eq. 

We note that the main quantities of interest within the 
PB theory can be expressed solely in terms of the cumu- 
lative density profile. This includes the PB potential field 
V^PB and the order parameters S^^. Using the definitions 
of these quantities f Section IIII|) . we obtain the following 
relations (which are valid for all ^) 



■0PB(y) = - J 



qPB _ ^ 



dye' 



npBiy')<iy' 
1 drtpB 



N dy 



(C7) 
(C8) 



3. Asymptotic behavior of S„ within PB theory 

a. Large Manning parameter ^ > Aap 

Consider the exact mean- field expression for Sf^^, Eq. 
(|40|l . The integrand in Eq. H4U|) is the product of an 
exponentially decaying factor with an inverse-squared 
sine-function, which has a series of divergencies at j/m — 
{mn — f)/ P for integer m and e = cot~^[(^ — I)//?]. For 
A — > cxD, we know from Appendix IC II that /? — > when 
^ > Aaf — 1 implying e — s- 0. In this limit, e may be 
expanded as 



P 



P' 



e-1 3(e-i) 



+ 0(/3^). 



(C9) 



The location of singularities, ym, tend to infinity with 
increasing A except for m = for which we have j/q = 
-e/P -1/(C - 1) using Eq. The quantity S^^ 

in Eq. H4()|l is therefore dominated by the lower bound 
of the integral (around y = 0) due to the exponentially- 
decaying integrand. To derive the asymptotic form of 
for large A, one can expand the integrand either 
around the lower limit of the integral y = or around 
the singular point j/o- Both procedures lead to the same 
scaling relation (|42|l for Sn in the strict limit of A — + cx3 
when ^ is close to the critical value = 1. But only 
the second procedure leads to a correct result when ^ 
increases beyond the critical value. This is because for 
large ^, the singularity at j/q ~ — 1/(^ — 1) approaches 
zero rendering the expansion around y — Q a poor ap- 
proximation. 

By expanding the integrand around y — Q (up to the 
first order in ?/), we obtain from Eq. H4UI) that 



PB 



n+l 



sm e Jo 



dye 



-(n+2i-2)y 



P^ + {^-l? 

^"+i(n + 2e-2)' 



(CIO) 



This relation reproduces Eq. H41|) in the text, which, as 
explained above, is valid for ^ close to = 1. 
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For larger values of we expand the integrand in 1401) 
around t/q — — e//3, which yields 



S. 



PB 



1 



(Cll) 



where r(a, b) is the incomplete Gamma function. This 
relation provides a quite good approximation for Sn for 
large A and in the whole range of ^ > = 1. In 
particular, when the limit A — s- cx) is taken, it yields the 
correct result for 5„(^, A — > oo) (see Eq. (jC18|) below). 



b. Small Manning parameter ^ < Aaf •' 

For ^ < Aaf — 1, the order parameters, Sn, vanish as 
A — > cx) indicating complete de-condensation of counte- 
rions. To demonstrate this, we use Eq. (|C8ll . which, for 
A » 1, can be written as 

= J^j^^y e""'^^PB(2/) + (C12) 

Since the cumulative density is bounded by the number 
of counterions, N , and tends to zero at any finite y for 
^ < Aaf ^ 1 (Appendix inSJ, we get S^^ ^ in this 
regime as A — > oo. 



4. PB solution in an unbounded system (A = oo) 

In the present study, we have assumed that the 
counterion-cylinder system is bounded laterally ensuring 
the normalization of density profile, ppb('?), to the total 
number of counterions, iV, even in the limit A oo. 
In a strictly unbounded system (with A = oo), the nor- 
malization property of density is not preserved, since a 
finite fraction of counterions escape to infinity. In this 
case, the PB equation H18() can be solved by relaxing the 
normalization condition l|21|l. Assuming the boundary 
conditions at the cylinder surface as tppni^) ~ ^ ^^'^ 
R[dij^Q{f = R)/df] = 2t one finds ^ 



21n- 



PB 



of counterions, au (Eq. (fZ^ ). i.e. 





e-1 C>e 



PB 



(C15) 



(compare with Eq. (|29|l ). The order parameters in the 
unbounded system, 8^^'°°, may be calculated using p^g. 
For ^ > = 1, we obtain 



jPB,oo 



1 



1 



r(o,- 



In the vicinity of the critical point (^ 1+), we get 

c 



(C16) 



^PB,oo 



(C) 



(C17) 



which exhibits a different exponent as compared with the 
quantity Sn^{(, A ^ oo) in Eq. This is again due 

to the difference in normalization factor, which enters 
in Sn through Eq. I|38|l (note the order in which the 
integration and the infinite-system limit are taken). In 
general, the order parameter S^^{(, A — > oo) is obtained 
by multiplying 5^^'°° (C) with the condensed fraction au , 
as 



5r"(C,A^oo)=aM^r'°°(C) 



(C18) 



APPENDIX D: HAMILTONIAN OF A PERIODIC 
CYLINDRICAL CELL MODEL (3D) 



2 In 



l + (e-l)ln. 



PB 



= 1, 

(C13) 
= 1. 



Also k2 = p^b(^) = for ^ < CP = 1 and = 
(C — 1)^/^^ otherwise. Hence using Eq. (|23|l . we have 
the density profile in a strictly unbounded system (for 
R<f<D) 



pTb(^1 









-2 r 


f 









1 -f (e - 1) In ■ 



e<i, 



which has the same form as given in Eqs. (|32|) and (|34|l . 

But now /5^(f) is normalized to the condensed fraction 



As stated in Section TV Bl the periodic boundary con- 
ditions used in the simulations in 3D lead to summation 
of Coulombic interactions (u3d(x) = l/|x|) over all peri- 
odic images. The resultant summation series are conver- 
gent for an electroneutral system and can be mapped to 
fast-converging series, which can be handled easily in the 
simulations 68, 69] . In what follows, we derive the con- 
vergent expressions for the configurational Hamiltonian 

The main simulation box (of height L and contain- 
ing N counterions) is replicated infinitely many times 
in z direction generating a series of M — > oo image 
boxes labeled with m = -M/2, . . . , -1, 0, 1, . . . , +M/2 
(with m = being the main box). The Hamilto- 
nian ||2Jl consists of three parts Hn — Ti-d + Hint + 
Tisoif, namely, counterion-counterion interaction, Tid, 
counterion-cylinder interaction, Hint, and the cylinder 
self-energy, Tisoif , that will be analyzed separately. (Here 
we use actual units and in the end, transform the results 
to the rescaled form.) 



Hint and HseU terms 



The counterion-cylinder interaction part per simula- 
tion box reads Hint / {M k^T) = S^=i u{ra), where using 
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(t(x) = asd{r — R), we have 

u{ra) = -q^B J V3b{x - XQ.)a(x)dx = 2^ In (j^^ + 

(Dl) 

with a running only over the counterions within the main 
box. The constant term is given by 

Co = -q^B J (x - xo)cr(x) dx, (D2) 

where xq belongs to the cylinder surface, cq may be 
written in terms of the cylinder self-energy, 



self «B 



^ J a(x)z;3D(x-x')a(x')dxdx'. (D3) 



Using the electroneutrality condition tL = qN (per box) , 
one can show that jSHscu/ {M kBT) = ~Nco/2. Thus, we 
have 

^ fWint + Wself) = 2^ ^ In f^) + Co, (D4) 



MkBT 



where Co = —PTCscu/M is a constant (see Eq. (0}), which 
diverges logarithmically with AI. This can be seen from 
the asymptotic behavior of the self-energy for large M 
(or large ML/R), i.e. 



dz dz'/{2ML) 



tHbML Jo 47r2 _ ^ 

where oq — In 2 — 1. This logarithmic divergency is can- 
celed by a similar divergent term coming from the inter- 
action between counterions (see below). 



2. Tici term and the Lekner-Sperb formulae 

The contribution from countcrionic interactions (per 
box and for large M) can be written as 



^ 9^ 

MfcsT 2M 



i-X!"3d(xj -Xj) 



(D6) 



2L 



X/3 



M 



where i and j run over all counterions (including periodic 
images), while a and (3 run over counterions in the main 

simulation box. We also have — 2X]m=^i ^^'^ 



■a — X^ 



A//2 I- -,-1/2 
2 , , 



E 

m=-M/2 



Pa/3 + (Ca/3 + n^) 



(D7) 



where pap = [{xa - xpY + {ya - ypfY^'^/L and = 
{za — zp)/L. Note that S^^j, in particular, represents the 
interaction between a counterion and its periodic images, 
which are lined up in z direction. This series diverges and 
may be written as 5°^ = 21n(Af/2) -I- 2Ce for M oo, 
where Ce = 0.577215 ... is the Euler's constant. In this 
limit, Sm is also divergent, but it may be split into a 
convergent and a divergent part as 



Sm — 



+ 5Ls(Pa/3,Ca/3) + 2(ln2-Ce), 



(D8) 



in which the convergent series S'ls(pq/3, Ca/s) can be ex- 
pressed in terms of special functions. 

Now inserting the above results for Sm and S^.j into 
Eq. (jD7p . we have for M ^ oo. 



q'h 



N 



J2 ^Ls(Pa/3,Ca/3)+e(Ce-ln2)+^7VlnM, 



MkBT 2L 

(D9) 

with a logarithmic divergent term (last term) from the 
one-dimcnsional periodicity of the system. Using Eqs. 
(|D5p and l|D9|) , it immediately follows that the divergen- 
cies in Eqs. (|D4|I and ljD9|l cancel each other when the 
electroneutrality condition is assumed. Thus we have the 
well-defined expression 



n-i ^ / 

I^-2e^ln(^- 

Q — 1 ^ 



q iB 

a#/3=l 



Shs{Pap,Ca0)+^Nho, 

(DIO) 

where ho = 1 + ln(i?/2L) + (Ce - In2)/7V; equivalently, 

Q=l ^ ' a^/3=l 

(Dll) 

in rescaled units, where /iq = 1 + ln(^^/2^Af) -I- (Ce — 
ln2)/iV. The above expression is used to obtain the in- 
ternal energy and the heat capacity of the system in our 
simulations (Section^). The term S'ls(Pq/3j Ca/3) niay be 
obtained from Eq. IID8|) using mathematical identities 
proposed by Lekner 68] and Sperb 69] . It may be writ- 
ten in the form of two formally identical series expansion 



I : - 21np„;3 + 4X;,'^=i iiro(27rmpc«/3) cos(27rmCa/3)- 



'LS 



II : y° 



(Pa/j) 



2m 



-Z(2m+l,l-Ca,3) 



Z(2m+ 1,1 + Ca/3) 
+ (P^/3 + S/3)-^/^ 



-*(! + Ca/3) -*(1- Ca/3) -C*. 

(D12) 

In the above relations, Kq(x) is the modified Bessel func- 
tion of the second kind, Z(n, x) — X^fc^o 1/(^ + 2^)" is the 
Hurwitz Zeta function, '^(x) is the Digamma function, 
and = 21n2 ~ 1.386 294. . .. 

The series in Eq. (jD12|l can be evaluated numer- 
ically up to the desired accuracy. Note that the 
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first series (Lekner sclienie) involves the Bessel func- 
tion Kq{x), which decays exponentially for large x (as 
~ exp{—x)/y/x) but diverges logarithmically for small x. 
It is therefore rapidly converging when the radial distance 
between two given particles, Pap, is sufficiently large. We 
use the following recipe to truncate series I: For paf3 > 3, 
we truncate after the third term, for 1/3 < pafs < 3, we 
include 2 + [S/pafs] terms in the sum (where [x] refers 
to the integer part of x), and for 1/4 < pap < 1/3, we 
sum at least 12 terms. This recipe ensures a relative 
truncation error of about \er\ ~ 10~^°. For small radial 
separation pap < 1/4 between two particles, series I be- 
comes inefficient and slow. We thus employ the second 
series expression (Sperb scheme). This series is rapidly 
converging for small pap provided that which enters 



in the argument of the Hurwitz Zeta function, is suffi- 
ciently small, namely for \(ap\ < 1/2 (note that in 
general we have \Cap\ < 1)- In fact due to the period- 
icity of the system, the energy expression (jD7|l remains 
invariant under the transformations (^ap — ^ 1 — Cap and 
Cap —Cap, and thus Cap can be always restricted to 
the range \Cap\ < 1/2. In this case, we use up to 8 terms 
in series II. The relative truncation error, \er\, varies for 
different Cap, e.g., for Cap — 0.4 and pap — 0.25, one 
has le^l - 10-^. The error substantially decreases for 
smaller Cap- The above truncation recipes yield accurate 
estimates for the interaction energies within the statisti- 
cal error-bars of the simulations. 
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